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Nonlinear  Scaling  Laws  for  Parametric  Receiving  Arrays 
Part  II  Numerical  Analysis 

by 

J.  W.  Kesner  and  F.  H.  Fenlon 


ABSTRACT 


This  report  outlines  the  procedures  used  in  a computer  program  for 
calculating  nonlinear  scaling  laws  for  parametric  receiving  arrays.  The 
basic  equations,  the  program  input  and  output,  the  program  listing,  and 
sample  computer  runs  are  included.  The  program  output  is  in  the  form  of 
curves  for  harmonic  amplitudes  as  a function  of  range,  the  "extra  dB  loss" 
as  a function  of  range,  and  parametric  receiving  array  beam  patterns  at 
given  ranges.  Both  detailed  analytical  solutions  and  approximate  solutions 
are  available  in  the  computer  program., 
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THEORY 


Part  I of  this  report  covers  the  tteoretical  analysis  that  forms 
the  basis  of  the  computer  program.  The  salient  equations  are  presented 

here  to  indicate  what  the  computer  program  does . 

The  spectral  amplitudes  »n  of  the  axial  pressure  field  are  given  by 

the  solutions  to  the  differential  equations 

* m 

J_(,  1 = -ML  - .(B-SWCN)t  - jD£  VmV 
t rr  v n j ■> 


where  , is  the  complex  spectral  amplitude  of  the  n harmonic  in 

n the  pressure  field 

R is  the  nonnalized  range  = range/Rayleigh  distance 

A = 0 for  plane  wave  cases 

= 1/R  for  spherical  wave  cases 

= R/ (1+R2)  for  mixed  (plane  and  spherical  wave)  cases 

B = 1/(1  +R2)  for  mixed  cases  with  phase 

= 0 for  plane,  spherical,  and  mixed-no  phase  cases 

C = n2(fs/f0)2(a0r0)ATTF 
£ /£  = source  frequency  to  pump  frequency  ratio 

„ ' ° = attenuation  coefficient  at  fQ  times  the  Rayleigh  distance 

0 0 2 2 
ATTF  = 1 + 0.69145*SALN/(l+n  *SWK  ) 

SALN  = water  salinity  in  parts  per  1000 

SWK  = sea  water  K 

D = (n/4)(fs/f0)oo 

0 = scaled  source  tone 

SWCN  = 0.69l45*SALN*DDT*Co*(koro)*n/(l+n2*SWK  ) 
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Co  = small  signal  speed  of  sound  (=1500  m/s  in  program) 

k r = wavenumber  times  Rayleigh  distance  at  f 
oo  u 

DDT  = 0.538*10'9*(1-6.54*10_4*PRES) 

PRES  = pressure  in  atmospheres . 


A differential  equation  solving  routine,  VODQ,  which  is  described  in 
Appendix  A,  is  used  to  solve  equation  (1)  for  4^.  Only  those  equations 
which  are  non-zero  are  included  in  the  solution  vector.  The  program  user 
specifies  the  number  (NHL)  of  lower  harmonics  (of  f ^ ) and  the  number  (NHH) 
of  higher  harmonics  (of  f ) that  he  wants  to  consider.  These  harmonics  are 
the  only  ones  that  are  used  in  the  differential  equation.  For  example  if 
NHL  = 3,  NHH  = 2,  and  f /f  = 10,  then  the  solution  vector  would  only 
contain 


V V V V V V ^lO*  Hll’  hi*  hi*  hi*  ^18'  ',19»  ¥20*  ^P  ^22’ 


By  definition  Y = 0. 


* 

and  v = 4*n  = complex  conjugate  of  4*^. 

By  specifying  NHL,  NHH,  and  f /f  the  infinite  sum  in  equation  (1)  becomes 
a finite  sum  over  at  most 


23’ 


NHL  + NHH*(2*NHL+1) 

different  positive  values  of  m.  For  harmonics  other  than  the  ones  deter- 
mined by  NHL,  NHH,  and  (fQ/ f g) v 4^  = 0. 

The  initial  values  of  v that  are  used  depend  on  the  type  of  problem 
that  is  being  solved.  If  the  user  simply  wants  to  see  the  harmonics 
generated  by  a single  frequency  then  the  initial  value  of  4^  = (0+^1)  and 
all  other  initial  4»  values  are  zero. 


If  the  user  wants  an  up-converstion  case  then 


’nfr  = (0+jX)  where  m " £o/£s 

and  = (O+.P) 

where  P = PSO  for  a plane  wave  case 

= PSO*EXP(-aor0/NFF.2)  otherwise. 

PSO  = signal  to  pump  amplitude  ratio. 

At  each  normalized  range  R the  incoming  signal  is  fed  into  the  differential 

equation  as 

= PSO*EXP(-ot  r *R/NFR2)*{sin(ARG)+-cos(ARG)} 

1 o o j 

where  ARC  = 2*(kQro) *R*(fs/fo)*sin2(0/2) 

and  0 = angle  of  intersection  between  signal  and  pump 

wave  normals. 

If  the  user  is  interested  in  a difference  frequency  case  then 

4NFR-1  = (°+j0’5)  311(1  ^NFR  = (°+j0’5) 

for  the  two  primary  frequencies  and  the  difference  frequency  amplitude 

will  be  generated  in  4^. 

The  extra  dB  loss  for  4^  is  computed  from  the  formulas 

EXDB  = - (20*logl0(4'1(R))  + 8.686*aor0*R}  for  plane  wave  cases 

EXDB  = - {20*login(4'(R))  + 8.686*o  r *(R-1)  + 20*log10(R)}  for  spherical 

iU  1 ° u wave  cases 

EXDB  = - (20*log10(4'1(R))  + 8.686*aQr0*R  + 20*log10(*nTR2)}  for  mixed  cases. 


APPROXIMATE  EXPRESSIONS 


Approximate  expressions  for  i^(R)  and  ’^■pR+i  (R>e)  were  generated 
in  part  I of  this  report.  These  are  also  included  in  the  computer 
program  as  an  alternative  to  solving  the  differential  equations.  The 
user  has  the  choice  of  two  expressions  for  i^(R).  One  uses  the  El 
function  and  the  other  involves  a definite  integral.  The  integral 
form  seems  to  give  better  results.  These  expressions  are  given  below. 
Approximate  expression  for  (R)  — El  form. 

^(R)  = EXP  (-aQro*R)  for  plane  wave  cases 

EXP  (-a  r *R) 

^(R)  = ° ° — 

Rp/^7/2)  2*EXP(4ctoro)*{El(2aoro)-El(2ctoro*RP)}2 

where  RP  = R for  spherical  wave  cases 
and  RP  = 1+R  for  mixed  cases. 

Approximate  expression  for  ^(R)  — integral  form. 


^(R)  = 


EXP  (-a  r *R) 
v o o . 


RP  / l+(a0/2)2*{£b  EXP(-2aoro*R)/RP  dR}2 


where  a=0,  b=R,  RP=1 
a=l,  b=R,  RP=R 
a = 0,  b = R,  RP  =v/l+R2 


for  plane  wave  cases 
for  spherical  wave  cases 
for  mixed  cases. 


Approximate  expression  for  ^j^pR+i(R»e) 

^NFR+l(R,e)  = 

{Ca0/2'»*EXP(-aNFR+1r0*R)/RP}*  fb  XP*^  (X)*EXP(aoro*X)*cos(KX)  dX  (2) 

3- 
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where  “NFRHro  * &W\ro 
K = 2(ksro)*sin2(6/2) 

ksro  - (WVo 
Vo  - (V>2/2 

k a = wavenumber  at  £ times  radius  of  circular  piston 
oo 

projector 

a - 0,  b = R,  RP  = 1,  XP  = 1 for  plane  wave  cases 

a = 1,  b = R,  RP=R,  XP=X  for  spherical  wave  cases 

a = 0,  b = R,  RP  = XP  =Vl+X^  for  mixed  cases 

0 = angle  of  intersection  between  signal  and  puirp  wave  normals. 
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PROGRAM  INPUT 


The  computer  program  may  be  run  from  a remote  terminal  or  by  cards 

as  a batch  job.  The  data  for  the  program  is  in  free  field  format  which 

means  that  the  data  does  not  have  to  be  in  any  particular  column. 

Input  variables  which  begin  with  the  letters  I,  J,  K,  L,  M,  N are 

integers  and  should  not  have  any  decimal  point.  All  other  variables 

must  have  a decimal  point.  For  lines  of  input  data  that  have  more  than 

one  number  the  data  values  must  be  separated  by  commas. 

Input  1 NOC  = number  of  separate  cases  being  run 

Input  2 UJP,  LPS,  NFR,  NHL,  NHH,  NAR,  NSG,  LAP,  NSC,  IBM 

LUP  = 1 for  one  frequency  harmonic  studies 

= 2 for  up-conversion  cases 

= 3 for  difference  frequency  cases 

LPS  = 0 for  plane  wave  cases 

= 1 for  spherical  wave  cases 

= 2 for  mixed  cases  with  phase 

= 3 for  mixed  cases  without  phase 

NFR  = fQ/f  = pump  frequency  to  signal  frequency  ratio 

NHL  = number  of  lower  harmonics  (of  f ) 

NHH  = number  of  higher  harmonics  (of  f ) 

NAR  = number  of  a r values 
o o 

NSG  = number  of  o values 
o 

LAP  = 1 if  the  approximate  expressions  are  to  be  used 
= 0 if  the  analytical  solution  is  to  be  used 
NSC  = 0 for  range  scaling  without  aQ 
= 1 for  range  scaling  with  oq 
IBM  = 1 for  beam  pattern  (only  if  LUP  =2  and  LAP  = 1) 

= 0 for  no  beam  pattern 


Input  3 
Input  4 
Input  5 
Input  6 


Input  7 
Input  8 


Input  9 


XMAX  = maximum  range  scale  value 

AR(L),  L = 1,  ...,  NAR  aQro  values 

SG(L),  L = 1 NSG  aQ  values 

SALN,  PRES,  TEMC,  ITA,  SWK,  AKDA 

SALN  = salinity  of  water  in  parts  per  1000 

PRES  = pressure  in  atmospheres 

TB1C  = temperature  in  degrees  centigrade 

(not  used  in  present  version  of  program) 

ITA  = 1 for  Marsh  $ Schulkin  form 

= 2 for  Russian  form 

(not  used  in  present  version  of  program 

SWK  = sea  water  K 

AKDA  = kQa  = wavenumber  at  fQ  times  piston  radius 
IAP  (input  only  if  LAP  = 1) 

LAP  = 1 for  El  approximate  form 

= 2 for  integral  approximate  form 
REM,  TMAX  (input  only  if  IK4  = 1) 

RBM  = normalized  range  R at  which  beam  pattern  is  taken 
TMAX  = maximum  angle  in  degrees  between  pump  and  signal 
normals.  The  beam  pattern  is  from  0°(on  axis) 
to  TMAX. 

PSO,  LPP  (input  only  if  LUP  = 2) 

PSO  = signal  to  pump  amplitude  ratio 
LPP  = 1 for  only  the  ,PNFR_1  curve  plotted 
= 2 for  only  the  ^TR+1  curve  plotted 
= 3 for  both  the  ii»NFR_1  and  'l'NFR+1  curves  plotted 


PROGRAM  OUTPUT 


The  program  output  is  in  the  form  of  curves  for  the  axial  harmonic 

amplitudes  (^(R),  02(R),  %r_i(r)»  %R+fR))  as  a £unction  o£  ran8e* 
the  "extra  dB  loss"  of  t^(R)  as  a function  of  range,  and  parametric 
receiving  array  beam  patterns  (R>6))  at  given  ranges.  A variety 

of  range  scales  are  used.  For  plane  wave  cases  the  range  scale  is  aQR 
or  R.  For  spherical  wave  cases  the  range  scale  is  cQln(R)  or  ln(R).  For 
mixed  cases  the  range  scale  is  oQ  arcsinh(R)  or  arcsinh(R) . The  "extra 
dB  loss"  curves  are  plotted  against  R.  The  ^npr_i(r)  ^ ^NFR+l^  curves 
are  plotted  against  R for  plane,  spherical,  and  mixed  cases. 

Many  cases  were  run  during  the  development  of  the  program  in  order 
to  test  the  validity  of  the  answers  that  the  computer  program  was  giving. 
Several  of  these  test  cases  are  included  in  the  figures  to  show  the  com- 
parisons between  different  assumptions,  such  as  plane  wave,  spherical  wave, 
and  mixed  wave  cases,  and  to  show  the  accuracy  of  the  approximate  expressions. 
Several  figures  are  included  to  show  the  versatility  of  the  computer  pro- 
gram since  it  will  handle  harmonic  studies,  up-conversion  cases,  and  beam 
pattern  studies.  All  the  figures  shown  were  also  chosen  with  the  idea 
of  demonstrating  how  the  harmonic  amplitudes  and  beam  patterns  change 
with  the  various  parameters,  such  as  aQro,  aQ,  and  fQ/fs. 

Figures  1 and  2 are  plane  wave  cases  that  show  i|^(R),  ^(X)*  and 
^3(R)  for  two  values  of  aQro  and  four  values  of  °0/a0r0* 

Figures  3 and  4 are  spherical  wave  cases  that  show  ij^(R),  ^(X)* 
and  extra  dB  loss  for  two  values  of  “0t0- 

Figures  5,  6,  and  7 are  mixed  wave  cases  that  show  ^(R),  ^(R)* 
and  extra  dB  loss  for  three  values  of  aQro  and  compare  the  curves  generated 
by  the  approximate  and  analytical  methods.  In  addition,  figure  7 covers  a 
range  of  five  0Q  values. 
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Figure  8 compares  the  spherical  and  mixed  wave  curves  for 
c A 4 with  similar  curves  previously  generated  by  Fenlon\ 

Figure  9 shows  the  similarity  of  the  two  up-converted  signals 

%FR-1^  and  *NFR+1^' 

Figure  10  compares  ^jyjpj^+^CR)  as  generated  with  the  mixed  phase  and 
mixed  no  phase  assumptions  and  with  the  mixed  approximate  integral 
expressions. 

Figure  11  shows  f°r  three  values  of  fQ/fs . 

Figure  12  shows  (R)  for  two  extreme  values  of  a r and  two 

6 NFR+1  o o 

extreme  values  of  a . 

o 

Figure  13  shows  up-converted  beam  patterns  for  three  values  of 
fQ/fs  at  two  different  ranges.  One  range  is  at  the  peak  of  the  ^Npp+i(R) 
curve  and  the  other  range  is  in  the  far-field  of  the  curve. 

Figure  14  shows  up-converted  beam  patterns  for  four  value  of  aQ 
at  two  different  ranges.  Some  of  these  patterns  are  near-field  patterns 
and  some  are  far-field  patterns  because  the  peak  of  the  curve 

changes  with  aQ. 


r m*m 
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CONSENTS  ON  THE  PROGRAM 

The  user  should  have  a good  idea  of  the  normalized  range  scale 
he  is  interested  in  so  that  he  can  input  an  appropriate  value  for  XMAX. 

It  was  found  that  using  NSC  = 0 (no  oq  scale)  was  often  more  convenient 
when  plotting  curves  for  a range  of  oQ  values.  For  up-conversion  cases 
using  the  approximate  expressions  it  is  important  to  use  a value  of  XMAX 
that  is  not  too  large  because  the  ^(R)  curve  is  precomputed  at  a finite 
number  of  points  over  the  range  up  to  XMAX  and  interpolation  is  used  for 
values  of  R other  than  those  prestored.  Bee a’  3 e of  the  steep  rise  and 
peak  of  the  (R)  curve  it  is  necessary  to  have  several  points  in  this 
region.  Otherwise,  a poor  representation  of  i^(R)  will  result  and  the 
resulting  curve  will  be  no  good.  Program  running  cost  was  the 

main  reason  for  prestoring  the  ^(R)  curve. 

One  suggested  change  to  the  program  is  to  include  a multiplying 
factor  Gm  inside  the  summation  of  equation  (1) , but  only  for  the  m = - °° 
to  m = 0 terms  in  the  mixed  with  phase  case.  The  suggested  Gm  term  is 

G m - 1/(1  + 2m/((l+R2)*(fo/fs))). 

Another  suggestion  for  making  the  program  more  complete  is  to 
replace  the  term  cos(KR)  in  equation  (2)  with  the  term  sin(KR)  + j cos(KR). 

If  the  product  KR  in  equation  (2)  is  large  enough  then  false  lobes 
can  appear  because  of  the  accuracy  limitations  in  numerical  integration 
of  a rapidly  varying  function.  These  false  lobes  have  shown  up  in  figure  13 
where  the  product  KR  is  very  large  for  some  of  the  cases. 
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As  an  added  feature,  the  computer  program  has  been  set  up  so  that 
difference  frequency  cases  can  also  e run.  Most  of  the  necessary  changes 
to  the  program  have  been  made,  but  nave  not  been  thoroughly  checked  out. 

If  anyone  would  like  to  use  the  difference  frequency  part  of  the  program 
it  is  suggested  that  he  contact  the  authors  first. 

For  more  details  of  exactly  what  the  computer  program  does,  the 
reader  is  referred  to  the  program  listing  which  is  included  in  this  report. 
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APPENDIX  A 


Included  in  this  appendix  are  descriptions  of  the  packaged  routines 
that  are  used  in  the  program.  The  routines  VODQ,  ROMBS,  and  DEI  are 
from  the  JPL  Fortran  V Subroutine  Directory,  Edition  No.  4,  October  1970. 
The  routines  TSCALE,  TSETUP,  and  TPLOT  are  from  the  Westinghouse  Research 
Laboratories  Fortran  library. 
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12.1.  initial  value  problem 


12.1.1.  NUMERICAL  SOLUTION,  S.P.  AND  D.P. 

12.1.1.1.  IDENTIFICATION 

NUMERICAL  SOLUTION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS. 

LANG  FILE  ELT/vERS  SIZE  ENTRY  NAMES 

F-V  LIU*JPLJ  VOIGO/JPL  5257 ( 8 ) =2735 ( 1C > VODQ, VODQl , VODQG 

SUBROUTINE  USED!  MNONE« 

COGNIZANT  PERSON:  F.  t.  KROGH,  SECTION  314,  JPL,  1969  JUNE  10 

12.1.1.2.  PURPOSE  - 

THIS  SUBROUTINE  COMPUTES  THE  NUMERICAL  SOLUTION  OF  THE  INITIAL 
VALUE  PROBLEM  FOR  A SYSTEM  OF  ONE  OR  MORE  ORDINARY  DIFFERENTIAL 
EQUATIONS. 


12.1.1.3.  REFERENCE 

1.  F.  T,  KROGH,  VODQ/SVDQ/OVDQ  - VARIABLE  ORDER  INTEGRATORS  FOR 
THE  NUMERICAL  SOLUTION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS, 
SECTION  314  SUBROUTINE  WRITE-UP,  JPL,  MAY  1969. 

2,  F.  T,  KROGH*  'ON  TESTING  A SUBROUTINE  FOR  THE  NUMERICAL 
INTEGRATION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS’,  JPL,  SECTION 
314*  TM  NO.  217,  MAY  1969. 


12.1.1.4.  PRECISION 

THIS  SUBROUTINE  IS  PRIMARILY  SINGLE  PRECISION,  USING  DOUBLE  PRECI- 
SION ARITHMETIC  ONLY  AT  A FEW  CRITICAL  POINTS,  TWO  OTHER  VERSIONS 
OF  THIS  SUBROUTINE  ARE  ALSO  AVAILABLE!  SVDQ,  WHICH  IS  ENTIRELY 
single  precision,  AND  DVDQ,  WHICH  IS  ENTIRELY  DOUBLE  precision. 
SEE  THEIR  WRITE-UPS  FOR  MINOR  DIFFERENCES  IN  USAGE. 


12.1.1.5.  REMARKS 

THE  ORDINARY  DIFFERENTIAL  EQUATIONS  MAY  BE  OF  ORDERS  1,2*3,  OR  4, 
AND  NEED  NOT  ALL  BE  OF  THE  SAME  ORDER.  IT  IS  SUGGESTED  THAT  IF  AN 
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ORDINARY  DIFFERENTIAL  EQUATION  IS  OF  ORDER  2,3»  OR  4,  IT  SHOULD 
BE  TREATED  AS  SUCH,  RATHER  THAN  BEING  REFORMULATED  AS  A SET  OF 
2,3,  OR  4 FIRST  ORDER  ORDINARY  DIFFERENTIAL  EQUATIONS,  THIS 
PERMITS  THE  USE  OF  INTEGRATION  FORMULAS  WHICH  HAVE  BETTER  NUMERI- 
CAL STABILITY. 


12.1,1.6.  METHOO 

THIS  SUBROUTINE  USES  LINEAR  MULTISTEP  PREDICTOR-CORRECTOR  FORMULAS 
OF  THE  ADA'S-MOULTON  TYPE  WITH  THE  APPROPRIATELY  GENERALIZED 
FORMULAS  FOR  ORDINARY  DIFFERENTIAL  EQUATIONS  OF  ORDERS  2,3,  AND 

4.  THE  SUBROUTINE  SELECTS  THE  ORDER  OR  9)  OF  THE 

FORMULAS  INDEPENDENTLY  FOR  EACH  ORDINARY  DIFFERENTIAL  EQUATION  OF 
THE  SYSTEM  AND  THEN  SELECTS  A STEP  SIZE  TO  BE  USED  FOR  ThE  ENTIRE 
SYSTEM.  THESE  SELECTIONS  OF  ORDER  AND  STEPSIZE  ARE  RECONSIDERED 
AND  POSSIBLY  CHANGED  AT  EACH  STEP.  THE  SELECTION  IS  BASED  UPON 
MAINTAINING  NUMERICAL  STABILITY  AND  MEETING  THE  USER’S  REQUESTED 
LOCAL  ACCURACY. 

the  subroutine  provides  special  returns  based  upon  either  the 

NUMBER  OF  STEPS,  THE  VALUE  OF  THE  INDEPENDENT  VARIABLE,  VALUES  OF 
THE  DEPENDENT  VARIABLES,  OR  VALUES  OF  AUXILIARY  FUNCTIONS.  ONlY 
THE  FIRST  TWO  OF  THESE  FEATURES  ARE  DESCRIBED  IN  THIS  WRITE-UP. 

SEE  REFERENCE  1 FOR  THE  OTHER  TWO  FEATURES  WHICH  INVOLVE  USE  OF 
THE  ENTRY  POINT  VODQG, 


12.1,1.7.  USAGE 


integer 

real 

real 

DOUBLE  PRECISION 


NEQ»KD(«N14) , IFLAG,MXSTEP,KSTEP»KEMAX, 
KQ(«N7H) 

EP(«N2«) ,HMINA,HMAXA,EMAX 
F(BN3tt) ,DT(lO»WN4tt) 

T,Y(*»N5W)  »H»DELT»TFINAL,YN(WN6tt> 


CALL  VODQ  (NE0,T*Y,F,KD,EP,IFLAG,H,HMINA*HMAXA,DELT,TFINAL» 
MXSTEf  KSTEP,KEMAX,EMAX,KQ,YN,DT) 

GC  TO  6 


4 CALL  VODQl  6 GO  TO  ( 10 , 10 , 3C , 4U , bQ ,60 , 70, 80 ) , I FLAG 


| 


01  m 


1?  [COMPUTE  F ( ) AS  A FUNCTION  OF  T AND  Y(  ).  IF 

IFLaG  = 2 The  value  of  t will  be  UNCHANGED  from  the  PRE- 
VIOUS RETURN.  THUS  THE  USER  CAN  EFFECT  SOME  REDUCTION  OF  EX- 
ECUTION TIME  BY  EVALUATING  ANY  SUBEXPRESSIONS  OF  F(  ) WHICH 
DEPEND  ON  T BUT  NOT  ON  Y(  ) ONLY  WHEN  I FLAG  = 1,  AND 


JPL  FORTRAN  V subroutine  DIRECTORY 


25  SEP  70 


12-4 


REUSING  THESE  VALUES  WITHOUT  RECOMPUTING  THEM  WHEN 
IFLAG  = 2.3 
GO  TO  4 


30  CAN  OUTPUT  POINT  HAS  BEEN  REACHED  UNDER  CONTROL  OF  DELT.3 
GO  TO  4 


40  C T HAS  REACHED  THE  VALUE  TFINAL.] 
GO  TO  . * . 


50  CAM  OUTPUT  POINT  HAS  BEEN  REACHED  UNDER  CONTROL  OF 
MXSTEP. J 
GO  TO  4 


60  CACCURACY  APPARENTLY  BEING  LIMITED  BY  ROUND-OFF  ERROR.  ' SUG- 
GEST THAT  EP ( ) BE  INCREASED. 3 
GO  TO  ... 

70  C ABS(H)  NEEDS  TO  BE  LESS  THAN  -HMlNA  TO  ACHIEVE  THE 
ACCURACY  SPECIFIED  BY  EP(  >.3 
GO  TO  ... 


«o  cerroneous  or  incompatible  setting  of  control  parameters. 3 


THE  DIMENSIONING  PARAMETERS  MUST  SATISFY: 

WNltf  /.GE.  1 IF  KD ( 1 ) >0 1 
N.GE.  NEQ  IF  KD(1)<0 

ttNCtt  /.GE.  1 IF  EP ( 1 ) >0 * 

\.GE.  IF  EP ( 1 ) <•)#  WHERE  «K«  IS  THE  SMALLEST  INTEGER 

FOR  WHICH  EP(ttKtt)  .GE.  0 

UN3H »«N4H#4N7«  .GE.  NEQ 


PN5«»ttN6tt  ,GE.  WKDTOTAL4 

WKDTOTAL4  /.EG.  KDUUNEQ  IF  KD(1)>0 

\.EQ.  SUM  OF  IABS(KD(D)  FOR  I = 1»...»NEQ  IF 
KD ( 1 ) <0 

BEFORE  CALLING  VODQ»  THE  USER  MUST  ASSIGN  INITIAL  VALUES  TO 


A4 
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T AND  Y ( ) Al'ID  SET  THE  PARAMETERS  nEQ,KD<  )*EP(  ),H,HMINa* 

hmaxa*uelt*tfinal,  and  mxstep.  the  user  may  change  the  values 

OF  EP { ) WHENEVER  IFLAG  = 1,  AND  MAY  CHANGE  THE  VALUES  OF 
HMINA»HMAXA,DElT,TFIMAL»  and  mxstep  at  any  time. 

THE  USER  MUST  COMPUTE  F(  ) WHENEVER  THE  SUBROUTINE  RETURNS  WITH 
IFLAG  = 1 OR  2. 

THE  SUBROUTINE  WILL  CHANGE  THE  VALUES  0^  T*Y(  >, IFLAG, H,KSTEP, 
KEMAX,EMAX*KG(  ) * YN(  ),  AND  DT ( ) DURING  THE  INTEGRATION,  ON 
RETURNS  WITH  IFLAG  r 3*4,  OR  5 THE  CURRENT  VALUE  OF  THE  SOLUTION 
WILL  BE  CONTAINED  IN  T»Y<  ),  AND  F(  >. 


THE  SUBROUTINE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS: 


NEQ  NUMBER  OF  DIFFERENTIAL  EQUATIONS  IN  THE  SYSTEM. 


INDEPENDENT  VARIABLE,  MUST  BE  SET  TO  ITS  INITIAL  VALUE 
BY  the  user,  thereafter  " T WILL  BE  CHANGED  BY  The 
SUBROUTINE. 


( Y( I ) , I = 1#«KDT01 ALtt)  DEPENDENT  VARIABLES.  FOR  SYSTEM  OF  FIRST 
ORDER  EQUATIONS  YU)  IS  THE  I**MTH#  DEPENDENT 

VARIABLE.  in  the  case  of  higher  order  equations*  the 

FIRST  LOCATIONS  OF  Y CONTAIN  THE  DEPENDENT  VARIABLE 
ASSOCIATED  WITH  THE  FIRST  EQUATION  AND  ITS  DERIVATIVES  UP 
TO  ORDER  ONE  less  THAN  the  ORDER  OF  THE  EQUATION  (WITH 
THE  LOWER  ORDER  DERIVATIVES  STORED  FIRST)*  THEN  FOLLOW 
THE  VARIABLES  ASSOCIATED  WITH  THE  SECOND  EQUATION,  ETC. 
FOR  EXAMPLE*  FOR  THE  SYSTEM: 

U* • = Fl (U,U*  *V,V* ,T> 

V**  = F2(U»U*»V,v*»T) 

THE  Y ( ) ARRAY  WOULD  BE  USED  AS  FOLLOWS: 

Y ( 1 ) = U 
Y ( 2 ) = U* 

Y < 5 > = V 
Y ( 4 ) = V* 


( F C I ) *I=1»NEG)  DERIVATIVE  VALUES.  IF  #DI«  DENOTES  THE  ORDER 

OF  THE  I**HTHW  EQUATION*  THEN  F C I ) IS  THE  8DI8**«TH8 
DERIVATIVE  OF  THE  I**«TH«  COMPONENT  WITH  RESPECT  TO  T. 
THE  USER  MUST  PROVIDE  THE  CODE  TO  COMPUTE  F(  )*  GIVEN  T 

i 


A5 
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AND  Y ( ) . 


KD(  I ORDER  OF  THE  DIFFERENTIAL  EQUATIONS  IN  SYSTEM.  DIFFEREN- 
TIAL EQUATIONS  WITH  ORDER  GREATER  THAN  4 CAN  ONLY  BE 
INTEGRATED  BY  BREAKING  THEM  into  LOWER  ORDER  equations* 
or  by  changing  certain  data  statements  in  the  integra- 
tor. IF  ALL  OF  THE  DIFFERENTIAL  EQUATIONS  ARE  OF  THE 
SAME  ORDER  THEN  IT  SUFFICES  TO  SET  KD<1>  EQUAL  TO  THAT 
ORDER  (KD(l)>0).  IF  DIFFERENTIAL  EQUATIONS  OF  MIXED 
ORDER  ARE  TO  BE  INTEGRATED*  KD  IS  A VECTOR  WITH  KD(1K0 
TO  INFORM  THE  INTEGRATOR  THAT  THIS  IS  THE  CASE.  THE 
ORDER  OF  THE  I**«TH«  EQUATION  IS  THEN  GIVEN  BY 
ABS(KD<I)),  (FOR  I>1  # KD ( I ) MAY  BE  EITHER  POSITIVE 
OR  NEGATIVE)  . 

EP(  ) parameter  used  to  control  the  LOCAL  ERROR.  if  ONE  wants 
THE  SAME  LOCAL  ERROR  BOUND  ON  ALL  COMPONENTS#  THEN 
EP ( 1 ) >0  AND  THE  ESTIMATE  OF  THE  LOCAL  ERROR  IN  EACH 

component  is  kept  less  than  EPm/i".  for  different 

ERROR  BOUNDS  ON  DIFFERENT  COMPONENTS#  LET  EP(1)<0  FOR 
I<WK«f  AND  LET  EP(MKtt)  .GE.  0#  THE  LOCAL  ERROR  CONTROL 
FOR  THE  I**«TH«  COMPONENT  IS  THEN  BASED  ON 
ABS(EP(I))  FOR  I<MKH  AND  ON  EP(«K«)  FOR  I .GE.  «Ktt. 
SUGGESTIONS  ON  HOW  TO  SELECT  EP(  ) ARE  GIVEN  IN 
REFERENCE  1. 


IFLAG  PARAMETER  USED  FOR  COMMUNICATION  BETWEEN  THE  INTEGRATOR 
AND  THE  USER.  THE  INTEGRATOR  SETS  IFLAG  AS  FOLLOWS: 

#< 

= i THE  VALUE  OF  Y C ) FOR  THE  CURRENT  STEP  HAS  BEEN  PREDIC- 
TED. THE  USER  SHOULD  COMPUTE  F(  ) AND  CALL  VODQ1.  IF 
A RELATIVE  ERROR  TEST  IS  DESIRED#  THE  NEW  VALUE(S)  OF  EP 

should  also  be  computed  here. 

=2  the  VALUE  OF  Y(  ) FOR  the  current  step  has  been  cor- 
rected. THE  USER  SHOULD  COMPUTE  F(  ) AND  CALL  VODQl. 

=3  AN  OUTPUT  POINT  HAS  BEEN  REACHED  (SEE  USAGE  OF  DELT  )* 

to  continue  the  integration  call  vodqi. 

=4  T=TFInAL.  IF  VODQl  IS  CALLED  WITH  T=TFlNAL  AND 
IFLAG=4,  IFLAG  IS  SET  EQUAL  TO  8.  IF  THE  VALUE  OF 
TFIIMAL  IS  CHANGED  THE  INTEGRATION  WILL  CONTINUE. 

=5  KSTEPzKSOUT  (SEE  THE  DESCRIPTION  OF  MXSTEP), 

=0  EMAX>C . 1 AND  IT  APPEARS  THAT  REDUCING  H WILL  NOT  HELP 
REDUCE  THE  GLOBAL  ERROR  BECAUSE  OF  ROUND-OFF  ERROR,  IF 
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THIS  OCCURS.  A LARGER  VALUE  OF  EP<  ) (OR  OF 
ABS  (EP  (KtNiAX ) ) IF  EP  IS  A VECTOR)  SHOULD  PROBABLY  BE 
USED.  IF  EP  IS  NOT  INCREASED,  TOO  SMALL  A STEPSIZE  IS 
LIABLE  TO  Bt  USED.  WE  HAVE  FOUND  THAI  REPLACING  EP(  ) 
WITH  32*EMAX*EP(  ) WORKS  REASONABLY  WELL.  (NOTE  THAT 
tMAX  IN  HOST  CASES  WILL  BE  ONLY  SLIGHTLY  LARGER  THAN 

o.l). 

=7  ABS1HXHMINA.  THIS  MAY  BE  THE  RESULT  OF  HALVING  H,  OF 
COMING  TO  THE  ENO  OF  THE  STARTING  PHASE  WITH 
ABS ( H ) <HM INA , OR  OF  THE  USER  INcREmSING  THE  VALUE  OF 
HMINA,  IF  ONE  WISHES  TO  CONTINUE  WITH  THE  CURRENT  VALUE 
OF  H,  SET  HMINA  .LE.  H AND  CALL  VODQl.  IF  THE  STEP- 
SIZE  HAS  JUST  BEEN  HALVED  (IN  WHICH  CASE 
EMAX>C • 1 ) , ONE  MAY  CONTINUE  WITH  THE  OLD  STEPSIZE  BY 
SIMPLY  CALLING  VODQl.  (SUCH  An  ACTION  IS  RISKY  WITHOUT 
A CAREFUL  ANALYSIS  OF  The  SITUATION.)  if  The  STEPSIZE 
HAS  I JOT  JUST  BEEN  HALVED,  SIMPLY  CALLING  VODQl  WILL 
RESULT  IN  THE  INTEGRATION  CONTINUING  WITH  THE  CURRENT 
VALUE  OF  H,  AND  A RETURN  TO  THE  USER  WITH  IFLAG=7  WILL 
OCCUR  AT  THE  END  OF  EVERY  STEP  UNTIL  ABS(H)  ,GE.  HMINA. 

=fl  ERROR  INDICATION,  THE  CONDITIONS  WHICH  MAY  CAUSE  IFLAG 
TO  BE  SET  EQUAL  TO  8 ARE  LISTED  BELOW, 


WHEN  CALLING  VODQl 

NEQ  .le,  n 
H*DELT  .LE.  Ol 

KD=e  OR  KD>H  (KD  A SCALAR)* 

KD(WIfl)=0  OR  ABS(KO(«I«))>4  FOR  SOME  » Itt »«Itt=l , 2 » . • . » 
NEQ  (KD  A VECTOR)* 


WHEN  CALLING  VODQl J 
IFLAG=4  AND  TrtFlMAL* 

EP=C  and  hmaxa*oi 

H*(TFlNAL-(  INITIAL  VALUE  OF  T))<i>  AND  THE  ESTIMATED 
ERROR  IN  EXTRAPOLATING  TO  TFlNAL  FROM  THE  INITIAL  POINT 
(USING  A FIRST  ORDER  METHOD)  IS  LARGER  THAN  PERMITTED  BY 
EPI 


WHEN  CALLING  VODQG  WITH  NG*0: 

NSTOP<0  OR  NSTOP>ABS<NG)  (SEE  USAGE  OF  THE  VODQG 
ENTRY  IN  THE  COMPLETE  WRJTE-UP) 


J 
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IF  VODQl  IS  CALLED  WITH  IFLAG  SET  EQUAL  TO  8,  THE  - 
PROGRAM  IS  STOPPED  AND  AN  ERROR  MESSAGE  IS  PRINTED. 


H THE  STEPSIZE,  H CAN  BE  EITHER  POSITIVE  OR  NEGATIVE. 

WHEN  SELECTING  THE  INITIAL  VAL^E  OF  H»  THE  USER  SHOULD 

remember  the  following: 

1.  THE  INTEGRATOR  IS  CAPABLE  OF  CHANGING  H QUITE 
RAPIDLY  AND  THUS  THE  INITIAL  CHOICE  IS  NOT  CRITICAL. 

2.  IF  IT  DOES  NOT  LEAD  TO  PROBLEMS  IN  COMPUTING  THE 
DERIVATIVES  (E.G.  BECAUSE  OF  OVERFLOW  OR  BECAUSE  OF 
TRYING  TO  COMPUTE  THE  SQUARE  ROOT  OF  A NEGATIVE  NUMBER )* 
IT  IS  BETTER  TO  CHOOSE  H MUCH  TOO  LARGE  THAN  MUCH  TOO 
SMALL. 


HMINA  MINIMUM  STEPSIZE  PERMITTED  (AFTER  THE  INTEGRATION  IS 

STARTED).  AFTER  THE  INTEGRATION  IS  STARTED,  AND  WHEN  H 
IS  HALVED,  AbS(H)  IS  compared  WITH  HMINA.  IF 
ABS ( H ) <HMINA  CONTROL  IS  RETURNED  TO  THE  USER  WITH 
IFLAG=7. 


HMAXA  MAXIMUM  STEPSIZE  PERMITTED.  THE  STEPSIZE  IS  NOT  DOUBLED 

if  doing  so  would  result  in  abS(H»>hmaxa. 


DELT  OUTPUT  INCREMENT • DELT  MUST  HAVE  THE  SAME  SIGN  AS  H. 

initially  tout  is  set  equal  to  the  INITIAL  value  of  t. 
whenever  t=tout  a return  is  made  to  the  user  with 

IFLAG=3.  (THUS,  IN  PARTICULAR*  A RETURN  IS  ALWAYS  MADE 
AT  THE  INITIAL  POINT),  WHENEVER  VODQl  IS  CALLED  WITH 
IFLAG=3,  TOUT  IS  REPLACED  WITH  T+DeLT.  IF  TOUT  DOES 
NOT  FALL  ON  An  INTEGRATION  STEP,  OUTPUT  VALUES  ARE 

obtained  by  interpolation  on  the  first  step  that 
(T-TOUT)*H>0.  INTERPOLATED  values  FOR  BOTH  Y AND  F 
ARE  COMPUTED. 


TFINAL  FINAL  VALUE  OF  T.  WHEN  T REACHES  TFINAL , CONTROL  IS 
returned  TO  The  USER  WITH  IFLAG=4,  OUTPUT  values  ARE 
OBTAINED  by  EXTRAPOLATION  IF  TFINAL  does  not  fall  on  an 
INTEGRATION  STEP.  IF  ONE  CHANGES  TFINAL  DURING  THE 
INTEGRATION  SO  THAT  (TFINAL-T>*H<0*  AN  IMMEDIATE 
INTERPOLATION  IS  PERFORMED  TO  OBTAIN  THE  OUTPUT  VALUES. 


MXSTEP  MAXIMUM  NUMBER  OF  STEPS  BETWEEN  OUTPUT  POINTS.  WHEN 
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VOOQ1  IS  CALLED  WITH  :*  .LE.  IFt.AG  .LE.  5,  KSOUT  IS 
SET  EOUAL  TO  KSTFP  4 MXSTEP.  A1  THE  END  OF  EACH  STEP, 
KSTtP  IS  INCREMENTED  AND  COMPARED  WITH  KSOUT.  IF 
KSTEP  .GE.  KSOUT  CONTROL  IS  RF TURNED  TO  THE  USER  WITH 
IFLAGrS.  THUS  IF  DELT  IS  SUFFICIENTLY  LARGE » CONTROL 
WILL  HE  RETURNED  TO  THE  USER  WITH  IFLAG=5  EVERY  MXSTEP 
STEPS. 


KSTEP  NUMBER  OF  INTEGRATION  STEPS  TAKEN. 


KEMAX  INDEX  OF  THE  COMPONENT  RESPONSIBLE  FOR  EMAX  (SEE  BELOW) 


EMAX  LARGEST  VALUE  IN  AMY  COMPONENT  OF  (ESTIMATED  ERROR)/«E« 

WHERE  WE«=ABS(EP)  FOR  THE  COMPONENT  UNDER  CONSIDERA- 
TION. ORDINARILY  THE  STEPSIZE  IS  HALVED  IF  EMAX>0.1. 
HOWEVERr  WITH  A RECENT  HISTORY  OF  ROUND-OFF  ERROR  LIMIT- 
ING THE  PRECISION,  VALUES  OF  EMAX  AS  LARGE  AS  1 ARE 

permitted. 


(KQ( I ) r 1=1, NtQ)  VECTOR  USED  TO  STORE  INTEGRATION  ORDERS.  KQ ( I > 
GIVES  THE  INTEGRATION  ORDER  USED  ON  THE  I**WTHW  EQUATION 


(YN( I ) , I=l,HKDTOTALK)  WORKING  SPACE.  (VECTOR  USED  TO  STORE  Y 
AT  THE  END  OF  EACH  INTEGRATION  STEP). 


((DT(I,J)»I=l»lG)rJ=l»NEQ)  WORKING  SPACE.  (DIFFERENCE  TABLES 
COMPUTED  BY  ThE  SUBROUTINE). 
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11.1.  ONE-OIMENSIONAI 


11.1.1.  QUADRATURE  , ONE -DIMENSION,  s.P. 

11.1.1.1.  IDEl  iT  IF  ICAT  ION 


QUADRATURE,  ONE-DIfOlSjONAL,  SINGLE  PRECISION 

LANG  PILE  Li  T/VERS  r,]7F  TNTRY  NAME'S 

E-V  LIB<<JPL<E  ROMPVJPl.  ITS!  ( 0 ) 1(10)  ROMPS, ROM?' 

SUBROUTINES  USED:  HNjnEM 


cognizant  persons: 


w,  R.  DIJNTON  AND  M.  UlFTHEl.M* 
JPL,  SECTION  314  » I PAr  SEPT  3n 


11.1.1.2.  PURPOSE 

OBTAIN  APPROXIMATE  EVALUATION  OC  A ONE -DIMENSIONAL  DEFINITE 
INTEGRAL.  By  NUmER rCAL  QUADRATURE, 

AMS  * THE  INTEGRAL  FROM  A TO  D OF  F(X)+OX 


11.1.1.3.  PI  FERENCrE J 


FOR  A COMPLE'S.  DESCRIPTION  OF  THIS  SUHROiiT  INF. . 

discussion  oc  a variety  or  test  cases,  see: 


INCLUDING  A 


WILEY  R.  BUNION,  MICHAEL 
QUADRATURE  SUBROUTINE  TOP 


UIETI'rt.M,  AND  I<ARFIJ  HAIGLEP, 
SINGLE  and  Ml.il  TIPI.E  INTI  OPAL 


’ROMBERG 
S’#  JPL 


w,  E MON, 
SUP  PC  JT  I fJs 
ME*«0P  UjBUM 


M. 

TO 

V 1 


DlETPEL"  • G . W 1 1 ! jr  , <MODJf  TED  PQMHr PG 
SUPPO  !T  GENERAL  SClLI'-'IflC  OQMpi  IT  JljG  ’ » 
3i4-:sn.  pp t | if  ,97.Jt 


OUADPATUPL:  A 

JPL  /MTEPUAL 


w.  DUNTON,  M o 
SUBROUTINES*  t 


D IE TH"|  t , *MC'P  j f AT’ONS  TO  T|)f 
TM  31 «>7,  l S[N-V4  J O', ,}. 


Jl’L  PGMIJFRG 


11.1.1.4.  METHOD 


THIS  SUBROUTINE  COMBINES  TECH  MlPi'rrr  r\)n  M 
STEP’  QUADRATURE  METHODS,  THE  ~G;  1,’i'Oi ' ' TMF 
SUBINTERVAL  CArOl.1  OF  THE  TOVA'  i ’ 

TO  APPROXIMATL  THE  INTEGRAL 


VOVA' 

cv;.,i 


1 1 


, v 


'ROMnERG* 
T f J i T 1 A L l 
i nvAi.  ca.hi 
SU'i  I NITP'/Al 


and  ’AOAPTIVE 
Y HICKS  A 
ANN  ATTEMPTS 
BY  USING  THREE 
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STAGES  OF-  ROMH-'IK,  QUADRA TUR(  . HUS  KEOU{|US  l VALUATION  OF  THE 
INTEGRAND  AT  F I VI  I 01 1 A | L V Sl’AGr  N POINT!.. 

IF  THIS  QUADRATURE  IS  REOai  DEI)  AS  SI  ICCEr'Sn  J|  (SEE  ‘ ERROR  CONTROL) 
THE  SUBROUTINE  Will.  ADO  T IK  VA(  UE  OBTAINED  TO  A RUNNING  SUM  AND 
PROCEED  TO  Tl  f AT  A NEW  DISJOINT  SUB INTERVAL  OF  THE  SAME  OR 
GREATER  LENGTH. 

IF  THIS  QUADRATURE  IS  IJOT  RI GARUFD  AS  SUCCESSFUL  AND  THE  CURRENT 
STEP  LENGTH  IS  GREATER  THAN  I IN 1 TO  VI  If  N "HE  SUBROUTINE  WILL 
REJECT  THE  RESULT  FOR  YMF  < UpRi;!K  SUN  INTERVAL  AND  TAKE  THE  LEFT 
HALF  OF  THE  SUblNTErwAL  AS  THE  I'EW  SUM I.MTf  RVAL  TO  HE  TREATED. 

IF  THE  CURRENT  STf  P LENGTH  T S Hft ) M Of’  SMALL. L R » THEN  THE 
SUBROUTINE  WILL  ACCEPT  TH!  CURpI  M"  !?rSUl-T  AND  pROCEFD  TO  THE  NEXT 
SUB  INTERVAL*  writing  a MESSAGE  ON  FORTRAN  UNIT  0 IDENTIFYING  THE 
SUBINTERVAL  0||  WHICH  THE  ACCURACY  TEST  WAS  NOT  SATISFIED. 


11,1,1.5.  error  control 

THE  QUADRATURE  IS  REGARDED  AS  SATISFACTORY  OVER  A PARTICULAR 
SUBINTERVAL  IF 


1.  THE  ESTIMATED  RELATIVE  E RP0R  0?  VIE  QUADRATURE  OVER  THAT 
SUB  1 1 ITER VAL  IS  AT  MOST  E UMAX  f OR 

2.  THE  ESTIMATED  ERROR  IN  THAT  Sljll  I NTS  RVAL  RELATIVE  TO  THE 
ACCUMULATED  VAJ  UF.  uF  T||E  INTEGE’AL  UP  TO  AND  INCLUDING  THAT 

SUB  interval  is  at  MOST  ermax,, 

3.  THE  ESTIMATED  aBSOINVL  ERROR  OVER  THAT  SUB  INTERVAL  IS  AT 
MOST  . 1«AbS(ERmAX) . 


11.1,1.6.  PARAMETER  CHECK  INC 


THF  SUBROUTINE.  WILL  TERMINATE  EXECUTION  WITH  A PRINTED  MESSAGE 
AND  A STANDARD  CXI.C  0 WALK  ••RACK  < PE  TURN  Gi  IT  TUI  GIVEN  PARAMETERS 
DO  NOT  SATISFY 


A<B  f On  A>B»  BUT  NOT  Ar B 
Oc<HMINcLE,HStAR.I.EcHmA)(  AND 


RELATIVE  - 0.<ERMAX<J  ,C»  BUT  NOT  FPMAX=0, 
ABSOLUTE  - ANY  MEGA  T I Vi  NUMBER  f BUT  NOT  =0. 


IE  A.GT.B  OR  IF  HSTAR.GE,  ONLY  A WARNING  MESSAGE 

IS  PRINTED. 


COPY  AVAILABLE  TO  DDC  DOES  NOT 
PERMIT  FULLY  LEGIBLE  PRODUCTION 
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11.1,1.7.  USAGE 

REAL  A*  U«  X.  FOFX.  Hf.TAR,  HMIN,  UMAX » ERMAX*  ANS 

integer  k»  key 

CaSSIGN  VALUES  TO  A»  D#  Mf,TAR.  MMlN,  UMAX*  ERMAX*  AND  Key] 

call  ROMBS ( a *D. X f FOFX.HSTAR.HMI N, UMAX, ERMAX* ANS,Kf KEY) 

10  [EVALUATE  THE  INTEGRAND  USING  THE  CURRENT  VALUE  OF  X AND 
STORE  THE  RESULT  IN  FOFX. 3 

CALL  ROM 2 

IF(K.EQ.I)  GO  TO  10 

C QUADRATURE  IS  COMPLETED.  RESULT  IS  IN  ANSI 


THE  SUBROUTINE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS: 


A,B  LIMITS  OF  INTEGRATION.  REQUIRF  A NOTr.H, 

X VARIABLE  SET  DY  THE  SUBROUTINE  FOR  INTEGRAND  EVALUATION 

IN  THE  USER’S  PROGRAM, 

FOFX  VALUE  OF  INTEGRAND  COMPUTED  BY  USER'S  PROGRAM  USING  THE 
ARGUMENT  X. 

H5TAR  SUGGESTED  INITIAL  STEP  SI7F.  THE  INITIAL  STEP  SIZE*  H* 
WILL  BE  SET  AT  11= , 0 J. * ( D-A  ) IF  UST AR  . OF. . ( B-A  ) /4  . OR  AT 
HrHSTAR  IF  HSTAR  ,L.T  . ( B~A  /\ . THE  FIRST  SUB  INTERVAL  WILL 
BE  OF  LENGTH  4.*H  AND  WILL  REQUIRE  FIVE  EVALUATIONS  OF 
THE  INTEGRAND  AT  X=A * A-M 1 » » . „ A*4*H.  SUGGEST 
HSTAR, GE, ( B - A ) / 4 , 


HMIN 

UMAX 


REQUIRE  hmin  ,le.  hstar  ,lf.  H»'AX. 

minimum  allowable  step  size,  reouiri 
hstar 

MAXIMUM  ALLOWABLE  STEP  SIZE,  REQUIRE 


0.  < HMIN  ,LE. 

hstar  ,le,  hmax. 


ERMAX 


TOLERANCE  ON  RELATIVE  OF  ABSOLUTE  ERROR,  SEE  DISCUSSION 
ABOVE  UNDEfI  'ERROR  CONTROL*.  REASONABLE  SETTINGS  FOR 
ERMAX  V/OUl  D BE  IN  the  RANGE  From  1,E.-.u  TO  l„E-7,  if 
GREATER  ACCURACY  IS  REQUIRED  , SEE  THE  i.'R I TF.-’UP * ON  ROMBD. 
IT  IS  REQUIRED  THAT  Oj-T .FRMX.l-T,. ) , F OR  THE  RELATIVE 
ERROR  TEST,  ONE  SHOULD  KNOW  THE  RANGE  OF  THE  ANSWER 


#.*s*rra'2'-*T 


JPL  FORTRAN  v SUBROUTINE  directory 


2S  SEP  70  11-5 


before  mi  i im's  aus oi. mi  irror. 

ANS  THE  riNAL  VAI.UF  0!  Till  IMTFGRAI  . AVAlUnLF.  WHEN  ROM? 

RETURNS  WITH  Kr;, 

K BRANCHING  PL  AO  SET  BY  TUf  SUBROUTINE  PO|<  IJSF.  IN  THE 

USER'S  PROGRAM , Krl  Mr. AUG  TIlE  USER  SHOULD  EVALUATE  THE 
INTEGRAND  AT  X,  STORE  THE  VAt.Ur  III  FOFX , AND  RE“ENTER 
R0M2.  Kr?  MEANS  THE  COMMUTATION  IS  fOMPLETED  AND  THE 
VALUE  IS  IN  AMS. 

KEY  FLAG  TO  CONTROL  PRINTING  OP  ERROR  MESSAGES,  PREVIOUSLY, 

ANY  VALUE  OF  KTY  NOT  = 7 WOULD  WRITE  A DIAGNOSTIC  MESSAGE 
WHEN  H BECAME. IT. HMIN,  THE  INPUT  VALUES  HAVE  BEEN 
CHANGED  SO  THAT  WHEN 

ACT  ION 

KEY  = S PRINT  INTERMEDIATE  T AND  Y VALlltSl  PRINT  THE  HMlN 
DIAGNOSTIC  IF  nETrCTEI). 

=0  PRIN1  jNTFRMf DlATl  T AND  y VALUES  1 DO  NOT  PRINT 
TNt  NMIN  DIAGNOSTIC. 

-7  DO  NOT  PRINT  Tllf  T AND  Y VAt  ULS  OR  HMIN 

diagnostic. 

= A|J Y OTHFR  VALUE  PRINT  TNF  HMIN  DIAGNOSTIC  IF 
DETECTED. 

THE  1 VALUFS  PRINTED  ARP  T (1  , 1 ) , AND  T<2.0).  THE 

PRINTEf  Y VALUES  ARE  THF  FUNCTIONAL  EVALUATIONS  Y ( 1 ) THRU 
Y ( 5 ) at  THE  .POINTS  X , X+N . . , . , X ! Oil . SEE  REFERENCES. 


U. l.l. a.  remarks: 

the  following  notes  may  re  applicable  for  hifficult  integrands. 

1.  if  in  doubt.  one  should  use  a small  value  or  hstar,  the  step 

SIZE  H CAN  DOUBLE  QUICKLY  AMD  THF.  USER  IS  PENALIZED  ONLY  A 

small  number  on  functional  (Valuations  wiuu  he  incrf.ases  nis 
chances  of  getting  an  accurate  APPROXIMATION  many  fold. 

2.  EJE  CAUTIOUS  WHEN  RELATIVE  1.RMAX  IS  .GT.  I"-o,  IF 

HSTAR.LT.  (B-AJ/L- , BUT  NOT  SMALL  ENOUGH . AND  THE  FUNCTION  IS 
OCCILATORY,  VERY  DIFFICULT,  ETC,.  ROMPS  CAM  RETURN  A WRONG 
ANSWER,  example:  FIX)=X*r.lU7,r>X*C0SX  ON  TUP  INTERVAL  (0.2PI), 

HSTAR  = A , 57  . Time  ANSW  R=™ . 2C9G724B . A RELATIVE  TOLERANCE  OF  .1 
GAVE  - . 25L-S * AND  A RELATIVE  ToLCRAMCP  OF  .fl  GAVE  A.lflOJ  DAD 
ANSWERS. 


i* 

mm 
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3.  ALSO  FOR  RELATIVE  F.RMAX  „ L T , 1 fl  * *-f, . IF  THE  RELATIVE  ERROR  IS 
ASKING  FAR  GREATER  THAN  6 SIGNIFICANT  DIGITS  ONI  Ir  PUSHING  THE 
ACCURACY  OF  THE  UNIvAf  1100„  ON  A HIGHLY  OCX  11  ATORY  # VERY 
DIFFICULT  PROBLEM,  ROMMS  MAY  1)1.  TAKING  THOU'., A! JOS  MORI'. 
FUNCTIONAL  EVALUATIONS  AND  1-jOf  lU'AUdNG  THE  At. CURACY  IT  DID  AT 

example:  same  f<x>  as  amove  when  gj\/i  n toi  erance  was 

10*4-6#  A N S = » <1 C ,r*  6 7 ?i  (\  0 WITH  A LOCO  E liNCT  lOi  IAI  iVAUATIOri'jJ  HUT 
WHCN  WHEN  F.F<MAX  = 10**-M»  ANS:: » 00067766  AND  10*267  FUNCTIONAL 

EVALUATIONS. 

4,  IF  ONE  WANTS  A ROUGH  APPROXIMATION  OF  THE  IMTRGRAL#  ME  CAN  SET 
HMIN=HSTAR=HMAX.  a Fixri)  svrp  INTEGRATION  of  the  function  will 
take  PLACE.  Bt  SURF-  THAT  l'|-.Y  = 7#  OR  MANY  DIAGNOSTICS  MAY  BE 
PRINTED. 


JPL  PORTRAN  V SUBROUTINE  DIRECTORY 
4.2.13.  EXPONENTIAL  INTEGRAL 
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4.2.13.1.  IDENTIFICATION 

exponential  integral 
lang  file  fi  T/UCDr 

F"V  Llb*JPL$  dei/jpl"  siZe  ENTRY  names 

subroutines  used:  dexp.dlog 


COGNIZANT  PERSON*  f w Mr  .n. 

SON.  E.W.  NO,  JPL,  SECTION  314,  1,70,  SEPT.  14 

4.2.13,2.  PURPOSE 
F0R  * 

(ExpX!T"/JuDf  *L  FR°M  T=-I^INITY  TO  T=X  OP 

:zvv^TrL.  r:  ei  c xf  Sxr^rsethe  m 

«*p l^rfSor  FR0M  T=Z  T0  op 


4.2.13.3.  METHOD 
4 • r?.  1 3 . 4 • ACCURACY 

FOUmnl  IcCURAC^ItAtISTICS^ERe'pounOj  VAC  1108  AN°  "* 


INTERVAL 
OF  x 


MAXIMUM  RMS 

relative  relative 


(-150,  -4)  2.20-16  5.1D-17 
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CF. 


(-4,  -1) 

7.50-17 

1.20-17 

(-1#  0) 

M. 70-18 

1 .40-18 

(?.  0.5) 

1 .60-16* 

3.L0-17* 

(?.ti.  6) 

5. 50-17 

1. CO-17 

(6.  12) 

1.6D-17 

3. 00"18 

<»2.  24) 

2. 90-17 

7.1D-18 

(24.  IOC) 

8,90-17 

1.9D-17 

.w.  ng»  comm. 

ACM  VOL. 

13.  «7. 

P P.  446-449. 


4,2.15.5*  USAGE 

double  precision  OEI,  X 

USE  THE  FOLLOWING  FUNCTION  IN  A FORTRAN  ARITHMETIC  STATEMENT: 
DEI(X) 

WHERE  UEI  IS  THE  COMPUTED  VALUE  OF  THE  EXPONENTIAL  INTEGRAL. 
4.2.15.6.  LIMITATIONS 

CKO) =- INFINITY  IS  APPROXIMATED  BY  -7.2D75  AMD  El (X  .GT, 
174.673)  IS  APPROXIMATED  BY  47. 2075. 

THESE  LIMITATIONS  WERE  ORIGINALLY  IMPOSED  FOR  THE  IBM/360  AND  HAVE 
NOT  BEEN  MODIFIED  FOR  THE  UNIVAC  11u8. 


4.2.15.7.  ERROR  EXITS 


NONE 


TSCALE 


Specifications : 

SUBROUTINE  TSCALE  (V,NPTS,VLOW,VHIGH, FIRST, IXORY) 
INTEGER  NPTS, IXORY 
REAL  V(NPTS) ,VLOW,VHIGH 
LOGICAL  FIRST 


Purpose: 

To  scan  the  array  V to  obtain  the  minimum  and  maximum  values, 
which  are  then  used  to  establish  a range  which  produces  a rational 
scale  for  the  axis. 


Usage : 

The  subroutine  should  be  called  as 


TSCALE(V , NPTS ,VLOW,VHIGH, FIRST , IXORY) 


where 


V = array  containing  data  to  be  scanned 

NPTS  = number  of  points  in  V to  be  scanned 

VLOW.HIGH  = receives  upper  and  lower  values  for  the  range  of  V 

FIRST  «*  .TRUE,  for  first  call  for  each  plot, 

.FALSE,  for  subsequent  calls  which  may  be  made 

for  a multiple  curve  plot.  Note  that  all  curves 

should  be  scanned  via  TSCALE  for  a multiple  curve 

plot . 

IXORY  = 0 for  abscissa,!  for  ordinate 
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TSCALE  (Cont'd) 


Note: 


The  abscissa  is  eight  grids  down  the  page,  and  the  ordinate 
is  5 grids  across.  TSCALE  need  not  be  used  if  the  range  of  values 
is  known. 
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Specifications : 

SUBROUTINE  TSETUP (XLOW ,XHIGH ,YLOW,YHIGH, IEORF, IGRID, LABELX, LABELY) 
INTEGER  IEORF , IGRID , LABELX ( 1) , LABELY (1) 

REAL  XLOW,XHlGH,YLOW,YHIGH 

Purpose : 

To  set  up  work  area  before  data  values  to  be  plotted  are  stored 
(via  TPLOT).  Labelling  information,  ranges  of  values  for  the 
abscissa  and  ordinate  variables,  and  options  for  grid  and  annotation 
are  input  parameters. 

Usage : 

The  subroutine  should  be  called  as 

TSETUP (XLOW, XHIGH.YLOW.YHIGH, IEORF, IGRID, LABELX , LABELY) 

where 

XLOW.XHIGH  = lower  and  upper  values  for  the  range  of 
(YLOW.YHIGH) 

the  abscissa  (ordinate)  variable 
IEORF  = option  to  set  the  format  for  annotating  the 
axes;  supply  1HE  for  E-format  (1PE10.3),  or 
1HF  for  F-format  (F10.3). 

IGRID  = option  to  select  grid;  supply  1HG  for  grid, 
lHb  (blank)  for  no  grid. 


T SETUP  (Cont'd) 

LABELX(LABELY)  = array  containing  the  information  to 

label  the  X(Y)  axis.  LABELX(l)  (LABELY (1) ) 
contains  the  number  of  characters,  and  the 
labelling  data  itself  is  stored  six  characters 
per  word  starting  in  LABELX(2)  (LABELY (2)), 


with  a maximum  of  30  (24)  characters. 


TP  LOT 


Specifications : 

SUBROUTINE  TPLOT(PLOT,ICHAR,NPTS ,XV,YV) 

INTEGER  ICHAR.NPTS 
REAL  XV(NPTS),  YV(NPTS) 

LOGICAL  PLOT 

Purpose : 

To  store  plot  data  and  initiate  plotting  via  PLOT  control. 

Usage : 

The  subroutine  should  bo  called  as 
TPLOT  (PLOT,ICHAR,NPTS,XV,YV) 

where 

PLOT  ■ . FALSE . if  data  values  are  to  be  stored  and 

plotting  deferred  (for  multiple  curves  on  one  plot) ; 
.TRUE,  if  data  values  are  to  be  stored  and  plotted. 
ICHAR  = single  character  used  as  plot  symbol  (e.g.,  1H*) 

NPTS  = number  of  data  points  to  be  plotted 
XV, YV  =■  arrays  containing  the  coordinate  values  to  be 

plotted,  i.e.,  (XV(I)) ,YV(I))  is  the  Ith  data  point. 

Restrictions : 

1.  If  more  than  one  data  point  occupies  the  same  plot  position, 
the  most  recent  value  processed  will  be  used. 

2.  Data  values  outside  range  specified  by  TSETUP  will  not  be 
plotted,  but  will  be  listed  below  the  plot. 
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Figure  1 

(pages  F2  to  F4  ) 


Plane  wave  case 

a r = 10" 4 

o o 

a / a r = 2,  5,  10,  100 
0 0 0 9 9 


It 


tv 


m « 


m m 
mm 


FI 


Figure  2 

(pages  F6  to  F8  ) 
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Figure  4 

(pages  F14to  F16) 
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Figure  5 

(pages  Fl8to  F22) 
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Figure  6 

(pages  F24to  F29) 
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Figure  7 

(pages  F31to  F38) 
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figure  8 
(page  F40) 
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Figure  9 

(pages  F42to  F43) 
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Figure  10 
(pages F4 5 to  F47) 
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Figure  11 
(pages  F49  to  F51  ) 
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(pagesF53  to  F56) 
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(pages  F58to  F63) 


Beam  patterns 
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Figure  14 
(pages  F65to  F66) 
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DOUBLE  PE  EC  I I DM  R>  Y<1£0>  > H» DELT* TFINAL*  T H • 120)  » TFMD»XMS'G» 


1 PL  A S T 

ti  ] HE  M 1 I OH  V *"6  0 1 ' I...I  >6  0 • » K D < 1 > » ER  1 > ? KO  <■  1 £ O')  j DT  < 1 0 j 1 £'  0 ) * 

1 c (i  £'0  1 . l 3 • 6u  > « HP  f'5)  * S-G 5 > » P 6 Cl  OP  < 8 0 > * YGQQD  <6U  > * YLAST  <60  > 


£ I Y M * 6U  ' 

COMMON  F'l.l  1 XP>  Xp*  P 3 1 ' PS2' XDB1 » I X»  IFL.ftG»  M j PS'Nl  ' F'SN8> 

1 RLO  < Uh'  > N • i.  \ '•  1 ' RAT*  THAX 

D I MEN ' 1011  OP  •'££?  5 1 • ' 0 •'  ££  ? 5)  * Pi  1 <82 ' 5 1 « F"S£  <22  »5>  > XDB 1 f'22  > 5 > * 
1 P'SHl  • 22 < 5 '•  ' P i M2  1 £‘2 » 5 > ? X T <’4 1 » 5 > » F*fi T <’41*5  ' 

C OMHDH  F l P 3'  T DPI' * P D f * T H « RLA ST  * RGOOD  ? YGOOD * YLRS  T * 


1 IVHM  Mi- 
ll 1 ME"  MS  I Or  I PL-'  •:  £2  :•  - PLY  * ££>  « F'LXT  <4 1 > * PLY  T -•  4 1 ) 

[I  I mem:  I DM  ICHR  • 5 )'•  LBXD  •£)  j LBXP  <4)  > LB  YD  4 ) « LBY1  <2>  * LBY2  <2>  « 
1 lB  Y f'l ' 3 > ’LB  T '3  < 3 '•  * L B Y 4 < -3  > j L B X T < 4 ) * L B Y T < 5 > 

COMMON  M« FPAT * ft  UP  0 > S I G 0 « F * LF'S* KS* NHL* MHL  £ 1 j NFP* NHL2* LUP* 


1 MSG?  I Pp < 1 SM * !•  SMI  ? KSN2*  AKOft*  IBM?  PBM 
10  00  POP  MAT  ■ '» 

1010  FORMAT  LUP=»  I.3*2X*5H  L PS=»  I3>  2X»  5H  NFR=*  I3*2X*5H  NHL=« 

1 I 3 * £’  • > 5H  NHH=*  I 3«  £X*  5H  NAP=  * I 3?  8 X * 5H  NSG='  1-3*  -''5H  LAP=?  1 8* 

£ £'X i 5H  NSC  = » I3'£X?5H  IBM=j  I-3-1 
1 0 £ 0 F □ P MAT  ' 5 H M 0 C = M ' 3 ) 

1 0 3 0 cOR HAT  •: 5H  P 3 0=? El  0.4?  2X ? 5h  LPP=  * 1 3> 

1 04  0 FORMAT  '.bH  XMAX  = » E 1 0. 4 •' 

1 0 5 0 PQFMAT  6H  A UP  0 - .«  5 < E 1 0 . 4 » 3X  '•  ) 

1 06  0 F 0 R M A T •'  6H  i 1 1?  l.i=  ' 5 '•  E 1 0 • 4 ? 3X)  ) 

1 070  FORMAT  ' t.H  * ALM=  * E 1 0. 4 j 3X  j 6H  F'RES='  E 1 0.  4 ? UX  ? 6H  TEMC=  ? F6.  1 ? -3‘X 
1 5H  I T H=*  1,3 1 37  - 5H  7 Ml.  = ? E 1 0 . 4 ? 3X > 6H  K i>A= * El  0.4) 

1030  FORMAT  '.5H  IAF'=M3) 

1 09  0 FORMA T ' 5H  P BM= » E 1 0 . 4 » 2X > 6H  TMAX= » F7 . £) 

£010  FORMAT  ■ 4H  E P= . E 1 0 . 4 - 3X. 7H  HMINR=» El  0. 4) 

8 08  0 FORMAT  «:88H  ERF  DP  )M  CONTROL  PARAMETERS":' 

8030  FORMAT  (4 OH  APPROX.  F.C'US.  HOT  VALID  FOR  DIFF.  FREQ.) 

£ 04 0 FORMAT  • 4 1 H APPROX.  ECU'S  . ONLY  VALID  FOR  MIXED  CASES) 

RE Fifi  10 fiO*  HOC 
PRINT  1 08 O' HOC 
DO  400  1=1 'HOC 

READ  1 0 m 0 * LUP * LP  3 * MFP * MHL » MHH * MAR * MSG >LAP* M SC  * I BM 

PRINT  1 01  O' LU^'  LPS * MFP*  NHl*  MHH*MAR* MSG' LAP' NSC> IBM 

L.i IR  = 1 FOR  HARMONIC  STUDIES' =2  FOR  UP-CGHVEPS  IOH»  =3  FOR  DIFF. 


FPEC. 

LP  3 = 0 FOR  PLANE  HAVE'  =1  FOR  SPHERICAL'  =2  •FDR  COMBINED  PH  AS: 
=3  FOR  COMBINED  NO  PHASE 

LAP= 1 IF  MANT  APPROXIMATE  EQUAT IONS' =0  OTHERWISE 
NSC  = 0 FOP  NO  SIGO  SCALE' =1  FDR  SI  GO  SCALE 
I BM= 1 FOP  BEAM  PATTERN' =0  OTHERWISE 
I BM= 1 ONLY  FOR  LUP =2 
INPUT  CAM  BE  ORGANIZE 
IF  <NSG  . ul . 5 -|  N'SG=5 
READ  1 000' XMAX 
PRINT  1040' XMAX 
R F h D 1 0 0 0 > (.  A Pi  (,  L ) ' L = 1 j 
R PINT  1 0 5 0 ' •:.  H R < L > > L = 1 j M A P ) 

READ  1 0 0 0 j • 3 G 1 L. ) »L=1  ? M S G > 1 


AND  LAP=1 
D BETTER 


, a^ussss 


L2 


f 

Li 


1530 
154  0 
1550 

156  0 C 

157  0 C 
1580  C 
1580  C 
1 6 0 0 C 
1 6 1 0 
166  0 

163  0 

1 64  0 


PRINT  1 06 0 , CSG  <L  > » L = 1 , NSG> 

PERU  1 000, SAI.N, PRES?  TEMC, ITft, SWF - HR OR 

PRINT  1 07  0,  SALN, PRES, TEMC, I TR» SW^, AR OR 

SALN=SALINI  TY  IN  PARTS  REP  1000,  =0.  FQP  FRESH  WATER 

PRES “PRESSURE  IN  ATMOSPHERES 

T E M C = T E M P E R R T IJ R E IN  DEGREES  CENTIGRADE 

ITA=1  FOR  MAPSH+SCHULRIN  FORM,  =6  FOR  RUSSIAN  FORM 

S Ml  = S E A lit  A T E P R » RR  0 A = K 0 ♦ A 

TDF  D=0. 68145^'SALN 

D0T=0.  538E-?^  Cl . 0-6. 54E-4^PPES •' 

IF  c.ITA  .EC.  1.)  TA=7.  267E-12M  0.  (1520.  0/  (TEMC+273.  0 •' 

IF  I T A . EC . £.)  TA=  c 1 . OE-3 ' 755 . 0)  (TEMC+273 . O')  *10.  0+*  < 934 . 0 


Ib5  0 

1 CT  SMC +£73 . O'  1 

166  0 

IF  clap  . NF.  1>  GO  TO 

167  0 

READ  1000, I AP 

168  0 

PRINT  1080, I AP 

169  0 C 

I AP= 1 COR  1ST  APPROX. 

1 7 0 0 

IF  IBM  .NE.  1 ' GO  TO 

1 7 1 0 

P E A D 1 0 0 0 « R P M , TMA X 

1 7£  0 

PRINT  1 090, PPM, TMAX 

17  30 

£05 

CONTINUE 

174  0 

IF  'LP S . GE.  £•'  GO  TO 

175  0 

PRINT  £040 

176  0 

GO  TO  999 

177  0 

£04 

IF  CL UR  . NE.  3'  GO  TO 

178  0 

PRINT  £030 

1790 

GO  TO  999 

1 8 0 0 

1 1 8 

CONTINUE 

1310 

IF  cl UP  .ME.  £)  GO  TO 

1 13 

FORM, =£  FOP  END  APPROX.  FORM 
£05 

£04 

1 18 

101 


1S£0 

READ  1000, PS 

0 , LRP 

18  3 0 

PRINT  1 03  0 , P 

SO,  LPP 

184  0 C 

LPF-1  FOR  ONLY  N-l 

185  0 101 

CONTINUE 

I86  0 

I CHR  C 1 ' = " 1 

187  0 

I CHR -:.£:•  =•••£  •' 

1880 

I CHR  C 3':.  - •"  3 ' 

189  0 

ICHRC4  • = • 4"' 

1 9 0 0. 

I CHR  5 > = ■ 5 / 

1 9 1 0 

LPXD  C 1 ;»  =1 

19£0 

LPXD  (£>  = R' 

1930 

IF  CLPS  .NE. 

0*  GO 

194  0 

IF  CMS C .EC. 

0-1  GO 

195  0 

LPXP C 1 > =6 

I960 

L P X P •"  £ :•  = •' I G IJ  ♦ R ' 

1970 

GO  TO  115 

198  0 £01 

LBXP  C 1 > =1 

1990 

L PXP  • £ :■  = R 

£ 0 0 0 

GO  TO  115 

£ 0 10  114 

IF  cLPS  .ME. 

1)  GO 

£020 

IF  CM SC  .EC. 

O'1  GO 

£ 0 3 o 

L PXP  < 1 1 = 1 0 

£04  0 

L PX  P £ ')  = ''  S I G 0 ♦ L •' 

£050 

LEXP  (3)  = 'NCR 

CURVE , =£  FOR  ONLY  N+l  CURVE, =3  FOR  DOTH 


TO  114 
TO  £01 


TO  116 
TO  £03 


d- 


L3 


£06  0 

GO  TO  115 

£07  0 

£ 0£ 

lb:>  p • i i =5 

£080 

LBXP  ' £ :•=  LN  < P > ' 

£080 

GO  10  115 

£ 1 0 n 

1 16 

IF  'NIC  .EC1.  0>  GO  TO 

£110 

LBXP ' 1 >=15 

£ 1 £ 0 

LBXP  ' £ :•  = ' 3 IG0»R'- 

£130 

LBXP  '■  3 > = PC'S  I NH 

£ 1 4 0 

L BXP  f 4 > = ' < P '• 

£ 1 5 0 

GO  TO  1 15 

£160 

£0  3 

LBXP  ‘ 1 :■  = 1 0 

£17  0 

LBXP  < £.)  = ' RRC  S I N " 

£ 1 8 0 

LBXP  * 3>  = 'H  <P>  • 

£ 1 9 0 

1 15 

CONTINUE 

££  0 0 

LBVD 1 1>  =13 

££  1 0 

L B V D ••  £ > = E X TRR 

it!  c'  c'  0 

LBVD  <3>  ='"DB  LOS"' 

££"30 

LBVD C 4 v = ' S' 

££4  0 

LBY1 - 1 > =4 

££5  0 

LBV  1 ••£>=' F"SI  1 

££6  0 

LBV£ < 1 > =4 

££70 

LBV£  <£  > = "PS  I £ '' 

££80 

L BYM • 1 > = 1 1 

££9  0 

LBVM  •:£.>  = ' BIFF.  ' 

£3  00 

LBVM  ■'  3>  = FREQ.  ' 

£310 

LBV  3 < 1 > =8 

£3£  0 

LBV3  <£.>  = "•  PS  I <CN-  '' 

£330 

L BV  3 f 3 > = ' 1 > x 

£34  0 

LBV4  <■  1 > =8 

£35  0 

LBV4  •',£.:•=  'PS  I - N+' 

£36  0 

LBV4-:'3>=  1>  " 

£ 3 7 0 

IF  < I EM  .HE.  1>  GO  TO 

£38  0 

LBXT cl) =13 

239  0 

LBXT  <£.>  =-"9NGLE  '- 

£4  0 0 

LBXT  3 '•  = "■  DEGREE'' 

£41  0 

LBXT  C 4 :•  = " S ■’ 

£4£  0 

LBVT  < \ ">  =£0 

£43  0 

LBVT  (£>  = • PS I < N + ■ 

£44  0 

LBVT = '1>  PRE ' 

£450 

L BY  T <4  > = 3 SURE  ' 

£460 

LBVT <5> = ' DB" 

£47  0 

1 £ 0 

CONTINUE 

£48  0 

DO  350  L=l»NflP 

£49  0 

9 OR  0=RP  a :> 

£5  00 

nn  351  M=1 ? NSG 

£5 1 0 

SIG0=SG  <rT> 

£5£0 

XMSG=XNflX/S I G 0 

£5  3 0 

IF  <MSC  .EC.  0>  XMSG= 

£54  0 

FRRT  = 1 . O -'NFR 

£55  0 

F£RP=FRRT ♦♦£  ♦H OP 0 

£56  0 

IF  •:'LUP  .NE.  £'-  GO  TO 

£57  0 

RK  OP  0=  (RfcOfl^E)  ••'£.  0 

£580 

FKR£=£ . O^FRfi T ♦RK OR 0 

. L4 


2590 

THET  = 0 . 0 

£600 

ST££=  - IN  (THET-  ;?.  0)  ) ♦♦£* 

£6  1 1 1 

1 02 

CONTINUE 

S' 62  0 

IF  (LHP  .NF.  1>  90  TO  119 

S6  2 0 

CHLL  HP P POX 

£64  0 

90  TO  351 

S65  O' 

119 

CONTINUE 

c'  *S  6 U 

NHL£‘=£'*NhL 

£'67  0 

NHL21=NHL£+1 

£630 

IF  CNF  R .LE.  NHL£>  60  TO  108 

£'69  0 

N=NHL+NHH4NHL21 

2700 

KSN=NHL£'+1 

£710 

DO  109  J=1>N 

£720 

NL J=NHL+1+ J 

£730 

Jf1N=MDD ' NL  J?NHL21> 

£74  0 

IF  < JNN  . EQ.  0>  JNN=NHL  £ 1 

£75  0 

Q JN=  >:  NL  J - 0 . 5>  / NHL 21 

£'7  6 U 

J DN=  I NT  CQ  JN » 

£ 7 7 0 

KF'=  JDN*NFP+  C JNN -NHL- 1 > 

£780 

KS  •:  J ) =i  R 

£790 

1 09 

CONTINUE 

£8  0 0 

60  TO  110 

£8 1 0 

1 08 

N=NHH»NFF'  + NHL 

2320 

KSN=NFP 

£33  0 

DO  111  J= 1 > N 

£34  0 

1 1 1 

K 3 C J > = J 

£35  0 

1 1 0 

CONTINUE 

£36  0 

1 SNl=k 'N-l 

£37  0 

KSN£=KSN+ 1 

£330 

NEQ=£»N 

£390 

- 

DO  105  K = 1 > N 

'£90  0 

V <K>  = 0 . 

£91  0 

U CK  :>  = 0. 

£92  0 

1 05 

CONTINUE 

£93  0 

IF  • LUF'  .NE.  £')  60  TO  103 

£94  0 

V f KSN  >=0.0 

£95  0 

ii.l  <K  SN  :*=1.0 

£96  0 

• 

60  TO  1 1£ 

£97  0 

1 0 3 

IF  CL  UP  .EO.  1>  1.1  Cl:- =1.0 

£980 

IF  CLUP  .NF.  3>  (30  TO  '112 

£990 

M CF S'N>  =0.5 

30  0 0 

l.l  (K SN 1 > =0.5 

3010  C 

ASSUME  PUMP  FREQS.  IN  NpR-l > NFP 

30£0 

112 

CONT INUE 

3030 

K [i  '•  15=1 

304  0 

NHN=  U 

3050 

IF  CLP 3 . NE.  0>  60  TO  106 

306  0 

c'  = 0 . 0 D 0 

3070 

PLO=0. 0 

308  0 

IF  CLUP  . E Q . £ > '-.I  C 1 = P S 0 

3090 

TF I NRL=XMS6 

3 1 0 0 

F P C 1 > = 1 . OE-3 

3110 

H-TFIMRLM . 0D-£ 

Ii  IFF. 


FREQ. 


in  1 


L5 


3120 

HMINA=TFINALM  . 0E--4 

313  0 

hmaxh=tfinal--2.  0 

3 1 4 0 

DELT  = TF INAL /2  0. L 

3 1 r i 

MX STEP =20 

31  f 0 

TFND=TF INAL-HMINA 

3 1 7 0 

GO  TO  107 

3130 

1 06 

IF  LLPS'  . GE.  2>  GO  TO  117 

3190 

P = 1 . U D 1 1 

32  0 U 

RLO= 1 . U 

32 1 0 

I P >:  L U P . E 0 . 2 > W < 1 > = P S 0 ♦ E X P < 

322  0 

TP  I NAL=DEXF'  ''X MSG'1 

3230 

Pp  C 1 - = 1 . OE-3 

3:40 

H=TF1NAL»1 . OD-2 

325  0 

HM I NA=TF I NAL^l . OE-4 

326  0 

H M A X A = T F I N A L • " 2 . 0 

32  7 0 

DELT=DEXP  c'XMSG/2  0.  0D0)  -1 . 0D0 

3230 

MXSTEP=20 

3290 

tfmd=tfinhl-hmina 

330  0 

GO  TO  1 07 

33  1 0 

1 1 '< 

R= 1 . OD-3 

332  0 

RLO=  0. 0 

333  0 

IF  CL U P .EC'.  2 > l.l 1 > = P S 0 ♦ E X P < 

334  0 

TP  INAL  = DE XP  CXM SG>  2 . 0 D 0 

335  0 

EP  C l > = 1 . OE-3 

336  0 

H = TFINALM . OD-2 

337  0 

HM I Nh=TF I NHL* 1 . OE-4 

333  0 

HNAXA=7 FINALS.  0 

3390 

DELT=0.  5D0*  CDEXP  YXMSG-'c'O.  00 0) 

-34  0 0 
34 1 0 
34c  0 

343  0 

344  0 

345  0 
34  *3  Li 
34  7 0 
343  0 
343  0 
I-:5  0 0 

35 1 0 

352  0 

353  0 

354  0 

355  0 
353  0 
3570 
353  0 
3530 
36  0 0 
36  1 0 
362  0 
3630 
364  0 


MXSTEF-20 

,.7  rM  jvftlUBSE  TO  000  DOES  NOT 

PERMIT  FULLY  LEGIBLE  PRODUCTION 

I YN  <Kl  ">  =0 
104  CONTINUE 
I X = 0 

CALL  VO  DU  <.NEC , R , y , F«  KOj  EF'*  I FLAG » H > HM I NA  * HMAXA » DEL  T»  TFI NAL  ? 
1 NX  3:  TEP  KSTEP  <•  KEN  AX  ? EMAX » KQ » YN « 0 T> 

GO  TO  6 
4 CALL  VO  DC  1 

6 GO  TO  '■  1 0 » 1 0 ? 3 0 ? 4 0 j 5 0 » 6 0 » 7 0 » S 0)  j I FLAG 
10  CONTINUE 

IF  '.i_Up  • NE.  2>  GO  TO  113 
C'EFR=P3  0*EXF'  ■: -F2AP»P> 

PRST=FFR£4R*ST2£ 

V < 1 ■'  = P E F R ♦ S I N C F P S'  T 

Y C2'>  =PEFR»COS  '•  FRS H 
113  CONTINUE 

CALL  DERI V <R * Y > 

GO  TO  4 

30  IF  R . GE.  TFND  > GO  TO  200 
CALL  OUTPUT  CR-  Y> 


L6 


365  0 

IF 

I X . GE . 

33) 

GO 

TO 

'3  0 0 

366  0 

IF 

<LP3  . EC1. 

0) 

GO 

TD 

4 

3670" 

IF 

•IPS  . L 3 ■ 

1 ) 

DELT  = 

9 ♦ >'  D E 

363  0 

IF 

a PS  . GE. 

3 

OELT  = 

1 F'  + 0 • 

3690 

i 3 

TO  4 

37 1?  0 

4 0 

GO 

TO  3 0 0 

371  0 

5 0 

CRL 

L OUTPUT  C 

P ' V 

> 

372  0 

GO 

TO  4 

373  0 

6 0 

EP  ' 

i)  = 33.  ♦EMRX*EP  ' 

1 ) 

374  0 

GO 

TO  4 

375  0 

70 

HMINH=HNINR/ 

1 0 . 

0 

•!  i ' 6 


! ? 9 U 
fi  fi 


F'F  I NT  £01  On  EP  Cl)  « HM I Nh 
TR1D=TFIHAL-HMINfl 
NHM=NHM+ 1 

IF  CNHM  . LE.  3)  GO  TO  4 
GO  TO  300 


xp  cxmsg-'So.  o r.i  o •'  - 1 . ooo) 

5 0 0 > ♦ cDE X P X M S' G •"  £ 0 . 00  0 ')  - 1 


33 1 0 

3 0 

PRINT  £030 

333  0 

GO  TO  999 

3330 

3 0 0 

CRLL  OUTPUT 'P*Y> 

334  0 

3 0 0 

CONTINUE 

335  0 

351 

CONTINUE 

336  0 

IF  a BN  .HE.  1>  .GO  TO  36  0 

'933'  0 

f RLL  T'FTI  ip  >:Ti.  0*  TMRXj  -50.  0*  0.  0»  1HF»  1HG*  LBXT  < LBV 

339  i.i 

DO  361  M= 1 * N SG 

39  0 0 

HO  363  K=  1 < 4 1 

39 1 0 

PLXT  ( K 1 =XT  M.) 

393  0 

■ 365 

PLYT  <K>  =PRT  (K*M> 

- 39  3 0 

IF  1 M . EP.  H;G)  GO  TO  363 

394  0 

CRLL  TPLOT < . FRL SE . > I CHP CM) >41> PLXT > P|_ YT ) 

395  0 

GO  TO  361 

996  0 

36  3 

C RL L TPL OT < . TRUE . * I CHP  CM ) >41* PLXT , PL  YT  ) 

397  0 

361 

CONTINUE 

3 93  0 

GO  TO  350 

399  0 

. 360 

IF  i'LUP  .HE.  1)  GO  TO  315 

4 0 0 0 

C RLL  T S E T U P < 0 . 0 > X M R X > 0 . 0 * 1 . 0*  1HE*  1 H G > L P X P*  LB ’i ' 1 ) 

4 0 1 0 

DO  3 o 1 N= 1 < NSG 

4 03  0 

DO  303  K= 1 > I X 

4 0 3 0 

PLX  1 K>  =XP  ■’  V > M) 

4 0 4 0 

303 

PLY  CK)  =P  S1  CK>M- 

4 0 5 0 

IF  C N . EQ.  NSG')  GO  TO  303 

4 06  0 

CRLL  TPLOT  (.FALSE. > ICHR CM) * IX> PLX*PLY> 

4 07  0 

GO  TO  301 

403  0 

" i 

CRLL  TPLOT  TRUE.  > I CHP  CM"-  , IX*  PLX*  PLY) 

409  0 

3 0 1 

CONTINUE 

4 1 0 0 

IF'  CLhP  .EC.  1)  GO  TO  317 

4110 

DO  304  M=  1 * NSG 

4130 

DO  305  K= 1 * I X 

4 1 3 0 

PLX  CK ) =XP ' L > M) 

414  0 

305 

PLYCK>=PS3CK'*M) 

415  0 

IF  • CM  . GT.  1 GO  10  9 06 

4160 

CRLL  TSChLE  CPLY? IX* VLQW* VHIGH* . TRUE,  j 1> 

__4170 

GO  TO  304 

L7 


CX3 


4180 

3 08. 

CALL  TSCALE  CPLY?  I X?  VLOW?  VH I GH?  .FALSE.  - 1 > 

4 1 9 0 

3 04 

CONTINUE 

420  0 

CALL  T SETUP  ' 0.  0?XMAX?  0.  0? VHIGH?  1 HE k 1HG?LPXP?L 

42 1 0 

DO  307  MM.NSG 

422  0 

DO  314  K=l? IX 

4230 

PLX  <■}■■■  ">  =XP  O'  . M.' 

4-4  0 

314 

PLY  <Y  > =PS2  <K . M:- 

4 £5  0* 

IF  0-1  .EC.  Ni'G'*  GO  TO  308 

42  60 

CALL  TPLOT  FALSE.  » ICHP  O'D  , IX?  PLX . PLY > 

4270 

b O TO  3 0 7 

4 c 80 

3 08 

CALL  TPLOT  < . TRUE  . > I CHP  <M)  ? I X ? PLX ? PLY:- 

4 £80 

307 

CONTINUE 

4 3 ill  1 1 

317 

DO  312  M=1 ? NS 6 

4 31  0 

no  310  K = 1.1 X 

432  0 

PLX  CK  ■ =XD<K»  N> 

4330 

3 1 0 

PLY  OO  =XDB1  CK ? M> 

4 34  0 

PLY  Cl>  = PLY  •:'£:■ 

4350 

IF  CM  . GT.  1 ■ GO  TO  311 

4 8 i.i 

CALL  T SCALE ' PLY* IX? VLOW? VHIGH? . TRUE.  . 1) 

437  0 

CALL  T SCALE  CRLX?  IX? PLOW?  RHIGH?  . TRUE.  ? 0.) 

4 33  0 

GO  TO  312 

4 380 

3 1 1 

CALL  T SC ALE CPLY? IX?  VLOW? VHIGH? . FALSE. ? 1> 

44  0 0 

CALL  TSCALE CPLX? IX? PLOW? RHIGH? . FALSE. ? 0> 

44  1 0 

312 

CONTINUE 

44  £0 

CALL  TSETUP • PLO? RHIGH? VLOW? VHIGH? 1 HE ? 1HG?  LBXD 

44  3 0 

DO  308  M= 1 ? NSG 

444  0 

HO  318  1=1.1-. 

445  0 

PLX  < K>  =XEi  CK  * M> 

448  0 

3 1 8 

PLY  O ) =XDB1  <K?M> 

447  0 

PLY  C 1 ) =F'LY  C 2 ') 

448  0 

I p CM  .EC.  NSG':*  GO  TO  313 

448  0 

* 

C ALL  T PLOT  C . FALSE . ? ICHP  CM > ? I X?  PLX ? PLY') 

450  0 

GO  TO  308 

4510 

313 

CALL  TPLOT  < . TRUE . ? I CHP  CM.)  ? IX?  PLX?  PLY) 

4520 

3 08 

CONTINUE 

45  3 0 

GO  TO  350 

454  0 

315 

IF  CL  UP  .ME.  '£.)  GO  TO  32  0 

455  0 

IF  CLPP  . EQ.  2.)  GO  TO  337 

458  0 

OO  331  M= 1 . NS  G 

457  0 

HO  332  K = 1 . I X 

458  0 

PLX  CK'>  =XH  cK  ? M '1 

458  0 

338 

PLY  CIO  -PSN1  CK  ? M :• 

48  0 0 

IF  CM  .GT.  D GO  TO  333 

48. 1 0 

CALL  T'SCAL  E CPLY ? I X ? VLOW . VH I GH  ? . TRUE . ? 1 > 

4820 

CALL  TSCALE  CF'LX?  IX?  PLOW?  RHIGH?  . TRUE.  ? 0) 

4830 

GO  TO  331 

484  0 

333 

C ALL  T SCALE  CPLY ? I X ? VLOW ? VH I GH . . FAL SE . . 1 > 

4850 

CALL  TSCALE  CPLX?  IX?  PLOW?  RHIGH?  . FALSE . ? O'.) 

4880 

331 

CONTINUE 

487  0 

CALL  TSETUP CPLO. RHIGH? 0. 0?  VHIGH? 1 HE ? 1HG.  LfcXB? 

488  0 

00  334  M= 1 ? NSG 

4880 

HO  335  K= 1 ? I X 

4700 

PLX CK) =xn  CK ? M> 

471  0 

385 

PLY  CK  1 =P3M1  *4  * M •< 

4 7c  0 

IF  CM  .69.  N 7 G ■ GO  TO  336 
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C HLL  1 PL  LI  T < . FPL  *E.  • K HP  -H-  , I X * PLX * PLY  ■ 

474  it 

GO  TO  334 

4 750 

3 ?i-. 

CRLL  TPLOT- .TRUE.  ? I LHP  CM  ? IX*  PLX*  PLY) 

4 7 T'  0 

334 

C GMT  I HUE 

4 7 7 0 ' 

Z'  » 

IF  ' LpP  .EC.  1 • GO  TO  350 

4 78  0 

DO  34  1 M = 1 * NS G 

4 /?  0 

DO  348  K = 1 * 1 X 

4900 

PLX  <r>  =XD  CK  - Mi- 

481  0 

348 

PLY  CK)  =PSH8  CK*  M:- 

4880 

IF  CM  . GT.  1 ■>  GO  TO  34  3 

4880 

CALL  TSCRLE  CpLY*  IX*  VLOW?  VHIGH*  . TRUE.  * 1) 

404  0 

ChLL  TSC RLE  C PLX * I X * PLOW * PH I GH * . TRUE . ? 0) 

4850 

GO  TO  341 

486  0 

843 

CRLL  TSCRLE - PLY* IX*  VLOW*  VHIGH* . FRLSE . * 1 ) 

48  7 0 

CRLL  TSCRLE  CPI  X*  IX*  PLOW*  P H I G H * . FALSE . ? O ) 

4880 

34  1 

CONTINUE 

4890 

CRLL  T SETUP - FLO* PH IGH* O. Q? VHIGH? 1HE* 1 HG? LEXD*  LBY4 

49  0 0 

DO  344  M= 1 * NSG 

4 810 

DO  345  1=1* IX 

498  0 

PLX  •: k =XD  CK?  M :• 

498  0 

845 

PLY  ) = PS N£  <K?M> 

494  0 

IF  CM  .EC.  NSG)  GO  TO  346 

4950 

CRLL  TPLOT  C . FALSE . * ICHR  CM » * I X*  PLX  * PLY) 

4960 

GO  TO  344 

4970 

846 

CRLL  TPLOT  C.  TRUE.  ? ICHR CM)  , IX?  PLX*  PLY) 

4980 

844 

CONTINUE 

499  0 

GO  TO  350 

5000 

380 

IF  CLUP  .ME.  3)  GO  TO  35  0 

5 0 1 0 

DO  384  M= 1 ? NSG 

5080 

DO  388  K= 1 ? IX 

5080 

PLX  OO  =XP  CK  * M> 

504  0 

388 

F'LY  CK)  =PS1  CK,  M - 

5050 

IF  CM  . GT.  1 :<  GO  TO  383 

506  0 

ChLL  T SCRLE  <PLY? IX?  VLOW?  VHIGH? . TRUE. * 1) 

507  0 

GO  TO  384 

508  0 

383 

CRLL  TSCRLE  (PLY?  IX?  VLOW?  VHIGH?  . FALSE.  ? 1 ■' 

5 t.i  9 0 

384 

CONTINUE 

5 1 0 0 

CALL  TSETUP  CO.  0?  XMRX?  <i.  0?  VHIGH?  1HE?  1HG?  LBXP?  LBYM,) 

5 1 1 0 

DO  381  M= i ? NSG 

5180 

DO  386  K=1*IX 

518  0 

PLX  CK.:-  =XP  (K?  M) 

514  0 

386 

PLY CK) =PS 1 CK * M) 

5 1 5 0 

IF  <M  .EC.  NSG)  GO  TO  385 

516  0 

CALL  TPLOT  < . CRLSE . ? ICHP <M) ? IX?  PLX? PLY) 

517  0 

GO  TO  381 

518  0 

385 

CRLL  TPLOT  c. TRUE. ? ICHR CM) ? IX*  PLX? PLY) 

5190 

381 

CONTINUE 

58  0 0 

350 

CONTINUE 

5810 

4 0 0 

CONTINUE 

5880 

MEND 

999 

END 
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DEL  T U=  TF  I NHL  -20.  000 
PEL  T 1 = TF I NRi_  1 . 0D4 
DEL  r2=TFl  UHl.  1 . PD3 
DEL  T 3=  T F I MHL  ' 1 . 0 D2 
DEL  T 4=  TF  I NHL  -'  1 . 0D1 
GO  TO  43 
p=i . emu 
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TF  I MftL=  DEXP  <XNSG  :■ 
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DEL 12= DEXP  '•  XM I G 1 . II  [1 2 * -1  . ODD 

1 . 0D2 > -1 
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THERE  hRF  OTHER  TEC HU I CUES  FDR  EVALUATING  THE  NESTED  INTEGRAL 
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IF  <’  LF'S  . E 0 . 1>  XX  = 1.0 
DO  100  1 R= 1 * 42 

XR=XX 

PLK  <’IR>  =XX 

p S I L K < I P •'  = r P S I < X P * H 0 R 0 * S I G 0 * L P S * I HP  ’> 

IF  '..PS ILK  '•  IR > . LE.  1.0E-2C"  PSIl.K  <IR>  =1  . OE-20 
IF  ':.IP  .LE.  10)  DEL  TI  = DELT  1 
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1510  IF  -IF'  .FT.  10  . PHD.  IP  . i_E . 20>  DELTI=DELT2 

1520  IF  (IP  .GT.  20  .HUD.  IP  .LE.  30:.  DELT  I=DELT3 

1530  IF  -IP  .GT.  30)  DELTI=PELT4 

154  0 IF  CLPS  .EC.  0>  DELX=I'ELT  I 

1550  IF  (LPS  . EQ.  1 > DELX=XX*DELTI 

1560>  IF  CLP"  . GE . £>  DEL X = C XX  + 0 . 5 > ♦ D E L T I 

1570  XX=XX+ OELX 

1530  100  CONTINUE 

159  0 3 0 IF  .-IBM  .HE.  1'  GO  TO  20 

1600  R=PBM 

1 6 1 0 T M X = T M Ft X ♦ S.  14 159 -l  8 0 . 0 

162  0 NTH -4 1 

1 6 3 0 GO  TO  d 1 

1640  20  NTH= 1 

1650  TMX=0. 0 

1660  21  DO  2 U 0 ITH=1 ? NTH 

1670  THET=  CITH-D  ♦TMX--40.  0 

1680  XT  CITH-.  M>  = THETM  8 0 . 0'"3.  14159 

1690  40  IX=IX+1 

1 7 0 0 IF  ( I p. M . N E . 1 > X D < I X » M > = R 

1710  IF  'LUP  . NE.  1>  60  TO  50 

172  0 P'S  I 1 =F  F'S  I ( P :>  fl  OR  0 ? S I GO?  LPS  > I HF') 

1 73  0 I F .'P'S  II  . IE . 1 . 0E-20  > PS  1 1 = 1 . OE-2  0 

174  0 PS1  <IX»M')=PSI  1 

175  0 IF  CL  PS  . GE.  ?.:•  GO  TO  52 

1760  IF  (LPS  .EQ.  1>  GO  TO  5,1 

1 7?  0 I F CNS  C . EQ . O ' XP  C I X > M - =P 

1 7 8 0 I F ('  N S C .EQ.  1 .>  X P < I X > M ) = S I G 0 ♦ P 
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.1810  51  IF  CNSC  . EQ.  O')  XP CIXj M> =DLD6  CP> 
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1930  IF  CLPS  . EQ.  0>  HP3M=1.0 

1940  IF  CLPS  . EQ.  1>  DPSM=R 

1350  IF  CLPS  . GE.  £>  DPSM=SQPT  Cl . 0+R^£) 
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20  00  A =0.0 

2010  IP  CLPS  .EQ.  D A=1.0 


Lll 


aoao 

a 03  0 
£04  0 
£050 
£060 
£070 


E-P 

H I T RP  = ''  t'-H  > ♦ 1 . UE-£ 
HM I M=  < E-R)  *1  . 0E-4 
HMRX=  < E-R ' ••'£.  0 
FPMRX= 1 . 0E-3 
KEY=0 


£08  0 

CRLL  POME? -'ft* 

E< 

X :>  FQFX 

. HS 

TRF:«  HM  IN.  HMftX*  EPMRX,  RNS»  KSTDF'j  ¥ 

£08  0 

6£  IF 

CLP 5 . EG. 

0 :■ 

DM  1 = 1 

. 0 

£ 1 0 0 

IF 

• LPS'  . EQ. 

1 ) 

PH  I =X 

£ 1 1 0 

IF 

'■  LF'S  . GE. 

£;- 

nni=:i 

OPT 

Cl . o +:■■♦♦£ ) 

£ 1 £ 0 

IF 

’■  X . G T . PL  ¥ ' 

1'))  GO 

TO 

1 04 

£ 1 8 0 

p:i 

C =P  S'  I L K - 1 

£ 1 4 0 

GO 

TO  103 

£15  0 

jl 

Cl 

rt 

PS  ILF  '•!  j 

£16  O 

PM=PLK  •:  1 - 

£170 

PS 1 C= 1 . 0E-£  0 

£ 1 8 0 

DO  101  I P=£ < 4£ 

.£  1 8 0 
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££5  0 
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1 01 

CONTINUE 
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£33  0 

IF  C ESTOP  . EQ.  1>  GO  TO  6£ 

£34  0 

GO  TQ  63 

£35  0 

61 

RN8=0.  0 

£36  0 

6 3 

F'PP=0QEF^RN7 

£37  0 

IF  Cl  EM  .ME.  l>  F'SNE  < IX?  M>  =RES  CF'PP) 

£ 3 8 O 

IF  <1  EM  .HE.  1>  GO  TO  60 

£38  0 

IF  CITH  .EQ.  1>  PRTO=RESCPPP':' 

£4  0 0 

PRT  C I TH  - M > =£ 0 . O^RLOG  1 0 CftBS  CPPR>  -"PRT O') 

£4 1 0 

IF  ''PRT  i I THk  M ) .LT.  -50.0)  FRT  CITH?  M) 

£4  £ U 

GO  TO  £ 0 0 

£4  3 0 

6 0 

CONTINUE 
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13 00  14  CONTINUE 
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TLN 1 =TLN1  -Y  ' £^L  1 > ♦K I ♦ S ••I T •-  •:  1 . 0+K  I ♦ ♦■'£  ♦ 3 :■ 

TL.  N£=-Y  ''£^L  I :•  ♦ 8RFF>k  I ♦♦£  ♦HTTP' 

TLN£  = TLN8+Y  < £♦'_  I - 1 > ♦K I ♦SWT  / • 1 . 0+L  I ♦♦£  ♦ i Wk  ♦♦£  ■' 
IF  ':.LC'S  .EC'.  0 > 60  TO  104 

IF  CLP '5  .EC.  1 > 60  TO  105 

IF  CLP 5 .EC.  £•'  60  TO  106 

TLN  1 =TL  N 1 -PR6 1 ♦Y  £*L  I - 1 > 

TLN6-TLN6-PF6 1 ♦Y  Cc>L I > 

60  TO  104 

1 06  Tl.Nl  = TLNl  -RF61  ♦Y  ' ,6  ♦ L I - 1 +P£H*Y  f.6*L  I > 

TLN£=  VLN6-RR6 1 ♦Y  «*L  I :«  -R£  1 I ♦Y  ''£♦  L 1 - 1 > 

60  TO  104 

1 05  TLN 1 =TLN1 -P I *Y C£*LI-1 > 

TL.N£=TI.  N8-R I ♦ Y «: £ ♦ L I 
104  CONTINUE 

F •:;£♦!_  J-1>=TLN  1 + TNI.  N1 
F •:  ;?♦!..  I - =TLH?+TNLN£ 

100  CONTINUE 
RLRS T=R 
RETURN 
EN  [i 


COPY  AVAILABLE  TO  DOG  DOES  NOT 
PERMIT  FULLY  LEGIBLE  PRODUCTION 
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1 1' 0 iUBPQU  T I ME  OUTPUT  1 P < V > 

1 1 0 DOUBLE  P P E C I S' I DM  P»  V 1 1 £ 0 > « PLN*  V M R 0 « X M 3 6 

1 2 0 COMMON  N»FPRT»R HP  n , 3 1 0 ij  • F > LP 1 < P I- , NHL * NHL £ 1 « NFP « NHL£ « LUP « 

1 ? 0 1 TT.G’  IRP?  PCM*  1 Ml  1 * P CNe' ? HP  UR  * I EM*  P'E'M 

1 4 0 D I ME N 3 ION  c C 1 £ 0 ) * P'S  •:  6 0 > 

1 5 0 C OMMON  B L 1 1 •••- X P - X D * P S 1*PC£>  X DEI  ? I X ? I FLRG * M » * S N 1 • P 3 N £ * 

16  0 1 RLDj X MR X * N S C - X T*PRT*TMR X 

1 7 0 D I MENS  I ON  XP  <’££ *5)  ? X D •-.££ * 5>  * PS  1 <22 ? 5'>  * PS£  C££ * 5)  > XDB 1 •:.££ * 5 ' * 

1 8 0 1 P S N 1 £ £ . 5 ) ? P S N £ < £ 2 * 5 ) * X T <41  * 5 ) * P R T ■ 4 1 » 5 » 

190  c'  0 9 U FOPMRT  1 7H  Ic:LRG=j  I U > 

£00  '3010  FOPMRT  1 3H  P = *EliJ.4«  ••  X * y H : I G U*P=  * E 1 11 . 4 * 3H*7H  LN  1 P 1 = • E 1 0 . 4 * 
£10  1 1 £H  "3  I Or  0 ♦ L N 1 p 1 = * E 1 0 . 4 * 3 X * 1 7 H 5 IG0»Rc'CS  I NH  1 P>  = * El  0.  4 :■ 

££ 0 3 0£ 0 FOPMRT  C3X*£H  N«5X*2H  V*  UX*£H  W»  1 IX*  4H  MRG?9X?9H  F.XTRR  DB* 

£30  1 5H  LOSS) 

£4  0 3030  fqrmRT  <£X*  13*  4 03 X > E 1 0. 4)  > 

£5  U 3 04  U FOPMRT  ••■'EH  3 I GO- * E 1 U.  4) 

£60  . IF  ax  .EC.  0>  PRINT  3040?  SI  GO 

£70  PRINT  £030? IFLRG 

£30  3P=P 


£9  0 

IF 

':‘NCC  .EC. 

1 •'  S P = C I G 0 ♦ S R 

3 0 0 C 

LPI 

=0  FOP  PL 

HNE  WAVE?  =1  FOR  SPHERICAL?  =£  FOP 

31  0 C 

- 3 FOP  COMBINED  NO  PHR 3 E 

320 

IF 

CP  . LT.  1 

. 0D-£0>  GO  TO  10 

3 3 0 

PLN 

=DLOG ( R ’ 

34  0 

IF 

( L P 3 .EC. 

0 • D B L = 8 . 6 8 6 ♦ R 0 R 0 ♦ P 

35  0 

IF 

CLF'S  .EC. 

1 > DBL=8 . 686^  CPLN+h  OP  i>  CP- 1 . CD  O')  ' 

360 

I F 

CL  PC  .GE. 

£ > D B L = 8 . 6 8 6 ♦ < R 0 R 0 ♦ P ■ + D L 0 G < D S 0 P T < 1 

370 

S0LP=PLN 

• 

33  0 

IF 

CNT.C  .EC. 

1 > S 0 L P = S I G 0 ♦ S 0 L P 

39  0 

GO 

TO  11 

400 

10  PLN 

=-9. 99D9 

410  i 0l_P=-9. 99E9 

4£0  11  CONTINUE 

4 3 0 S 0 G R = D L 0 G < R + D 3 Q R T < p ♦ ♦ £ + 1 . 0 D 0 > 

44  0 IF  •: N S C .EC!.  1 > S 0 G R = S I G 0 ♦ S 0 6 R 

450  IF  'IFLRG  .EC.  5>  GD  TO  16 

460  IX=IX+1 

47 U IF  L P C .EC.  U>  XP  > I X * M ' =SP 

480  IF  <LPS  .EC.  1>  XP  < I X » M>  =S  OLR 

490  IF  <LF'C  . GE’ . £ ■'  XP  '■  I X?  M>  = i UGR 

500  VD<IXjM')=R 

510  16  CONTINUE 
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53  0 C 

IF  1 Mi-ij  • EU.  1 1 PR  I NT  3 U 1 ij  * R , i'R?  Rt_N,  U L R , CiJGR' 

530  C 

PRINT  3030 

54  Ci 

DO  301  K I = 1 , N 

550 

K-R  S O-  I > 

56  0 

y 21=2^1-1 

57  0 

K2=2*KI 

53  0 

V M 9'3 = D . i. ! R 1 < Y K 2 1 1 ♦ ♦ 3 Y ' . K 2 '•  ♦ ♦ 3 ■' 

59  0 

IF  <YMAG  . LT.  1.0  0-3  0 .DP.  P . LT.  1 . OH -3  0"'  30 

6 Ci  0 

E D E L = - 3 0 . 0 ♦ D L 0 G 1 0 C Y M H 6 > - D E L 

6 1 0 

'30  TO  1 3 

6 c!  0 

13 

EDEL=  9 . 99 E 9 

63  0 

13 

CONTINUE 

64  0 

IF  <1  FLAG  .EC!.  5 •'  '30  TO  14 

65  0 

Ip  '■  K I . NE.  1'  GO  TO  15 

66  0 

c'  5 1 < ! X , M •'  = Y U **  '3 

67  0 

YOPl  •'  I X»  M>  =Et'EL 

63  0 

GO  TO  14 

69  0 

15 

IF  <PI  .NE.  3)  GO  TO  17 

7 0 0 

p S 3 1 I X » M > — Y M A G 

710 

GO  TO  14 

730 

17 

IF  <K I . NE.  K SN 1 :•  60  TO  18 

73  0 

F'iNl  <IX,  M)  = YMhG 

74  0 

GO  TO  14 

75  0 

18 

IF  (KI  .NE.  KSN2)  GO  TO  14 

76  0 

IX,  M>  = YMAG 

77  0 

GO  TO  14 

73  0 

14 

CONTINUE 

79  0 

IF  <KI  . GT.  3>  '30  TO  301 

3 0 0 C 

PRINT  3 03  U » k , Y '.pir’l  •'  , Y ' K c!  •'  , YMRG , ED  EL 

31  0 

301 

CONTINUE 

S3  0 

RETURN 

33  0 
S'?  END" 

END 

